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Preface 


( 

This  report  contains  the  results  of  my  efforts  to  develop  and  de- 
scribe a relatively  fast  and  reasonably  accurate  computer  code  designed 
to  simulate  time-dependent,  high-pressure  shock  propagation  through 
a gas  or  solid.  The  most  important  result  of  my  work  is  the  one-dimen- 
sional hydrodynamic  computer  code  HUFF  that  I developed.  I have  tried 
to  present  a complete  description  of  the  derivation  of  the  model,  the 
assumptions  made,  and  the  steps  necessary  to  define  and  solve  problems 
using  the  code  here  at  AFIT.  I made  an  extensive  effort  to  provide  a 
simple  and  flexible  program  so  that  HUFF  could  be  used  in  whole  or  in 
part  in  solving  related  problems  in  the  future. 

There  are  many  people  to  whom  I am  indebted  and  would  like  to 
( thank.  First,  1 would  like  to  thank  Captain  Clarence  Anderson  for  pro- 

viding me  with  numerous  references  and  extensive  background  to  begin 
ny  work.  I would  Also  like  to  thank  #y  advisor,  Major  George  H.  Nickel, 
for  his  help,  guidance,  and  encouragement  during  the  research  period. 
Last,  but  most  important,  I would  like  to  thank  my  wife  Martha  for  her 
patience,  love,  and  understanding  throughout  the  period  of  this  work, 

David  J.  Peters 
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Abstract 


IIUFF  is  a one-dimensional  Lagrangian  hydrodynamics  computer  code 
developed  from  the  basic  principles  of  mass,  momentum,  and  energy  con- 
servation for  strong  shock  propagation  in  a solid  or  gas . Two  equa- 
tions of  state  n re  usod  - the  adiabatic  ideal  gas  law  with  a variable 
gamma  and  the  Gruneisen  solid  equation  of  state  with  a constant 
Gruneisen  ratio.  The  ftichtmyer  and  Merton  difference  equations  for 
strong  shocks  are  used  on  a spatial  mean  composed  of  up  to  100  cells. 
Results  for  two  problems  are  presented  which  show  the  usefulness  and 
limitations  of  the  code  and  also  serve  as  sample  problems.  The  results 
of  a one  kiloton  nuclear  explosion  are  compared  to  the  Nuclear  Blast 
Standard  (1KT) . The  results  were  within  13  percent  for  shock  over- 
pressure and  overdensity,  5 percent  for  shock  material  velocity,  and 
2 percent  for  shock  position  over  a range  of  20  meters  to  2 kilometers 
from  the  burst  point.  The  larger  deviations  occurred  at  early  times 
being  attributed  to  an  absence  of  radiation  transport  calculations  in 
the  code.  The  second  problem,  a megabar  compression  of  uranium,  shows 
agreement  within  two  percent  for  all  parameters  (peak  shock  pressure, 
density,  material  velocity  and  shock  velocity)  when  compared  with  the 
Rankine-Hugoniot  compression  curves . The  equation  of  state  for  a solid 
was  limited  to  calculations  below  100  megabars  due  to  its  simplicity 
and  constant  value  for  the  Gruneisen  ratio.  A complete  users  guide 
and  program  listing  are  also  provided. 
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HUPF , A ONE -DIMENSIONAL 


HYDRODYNAMICS  CODE  FOR 
STRONG  SHOCKS 

1 • Introduction 

One  of  the  primary  effects  of  low  altitude  nuclear  weapon  deto- 
nations is  the  development  and  propagation  of  strong  shock  waves . 

Large  computer  codes  ha>re  been  developed  to  accurately  describe  the 
development  and  propagation  of  shock  waves  in  air  and  in  solids  (Ref 
10 tl).  These  computer  codos  generally  require  large  computer  storage 
and  long  computation  times. 

In  1975,  Clarence  L.  Anderson  developed  a computer  program  based 
on  conservation  of  mass,  momentum,  and  energy  to  calculate  hydrodynamic 
shock  propagation  in  solids  and  in  air.  Anderson's  program  used  a 
course  spatial  mesh  (100  cells)  with  a rezoning  capability  and  some 
simplifying  assumptions  to  develop  a less  accurate  but  more  simple 
operating  code  (Ref  1:1).  The  results  of  Anderson's  work  was  encouraging 
and  lead  to  the  development  of  HUFF. 

The  primary  purpose  of  this  thesis  was  to  develop  a simple,  stable 
school  use  hydrodynamic  computer  code  based  on  conservation  principles 
and  including  appropriate  simplifying  approximations  providing  reason- 
ably accurate  results  with  relatively  short  computational  ti.me  require- 
ments . 

This  report  presents  the  conservation  principles  on  which  the  nora- 
( puter  code  is  based  and  the  application  of  these  principles  to  the 
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computer.  Development  of  the  computer  code  is  discussed  and  the  appli- 
cable derivations  are  presented.  Hydrodynamic  shock  propagation  in 
air  and  in  solids  is  developed  and  sample  problems  are  presented. 

The  development  of  the  hydrodynamic  shock  code  is  presented  in 
Chapter  II.  Chapter  III  discusses  the  initial  conditions  and  boundary 
conditions  applicable  to  the  HUFF  program  and  their  coupling  with  the 
HYDRO  subroutine.  Chapter  IV  discusses  the  results  of  two  sample 
problems  and  Chapter  V contains  the  conclusions  of  this  study. 


II . Hydrodynamic  Shock  Propagation 


This  chapter  presents  the  basic  differential  equations  that  govern 
the  propagation  of  strong  shocks  in  both  gases  and  solids.  Various 
simplifying  assumptions  are  made  to  reduce  the  equations  to  a one-dimen- 
sional fora.  The  Lagrangian  coordinate  system  and  subscript  notation 
used  in  HUFF  are  presented  followed  by  the  finite  difference  scheme 
used  to  approximate  the  governing  equations.  The  stability  conditions 
that  constrain  the  size  of  the  time  step  and  the  variable  gamma  modifi- 
cations co.  ~ ■ jte  the  chapter. 

Basic  Governing  Equations 

The  basic  equations  governing  strong  shock  propagation  are  derived 
from  the  physical  laws  requiring  that  mass,  momentum,  and  total  energy 
be  conserved.  For  propagation  of  strong  shocks  in 'solids  the  hydro- 
static pressures  are  very  large  and  the  solid  can  be  treated  as  a fluid. 
Thus  the  differential  equations  expressing  mass,  momentum,  and  total 
energy  conservation  apply  equally  to  both  gases  and  solids  (Ref  5:7). 

(u*7)p  + pV*u  ■ o (1) 

(u»7)u  + i VP  »*  o (2) 

P 

(u*7)E  *;•  — 7*(Fu)  ■ 0 (3) 

P 


The  equations  are: 
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where  p is  density,  u is  velocity,  b is  time,  p is  pressure 
and  E is  energy. 

The  equation  of  state  for  the  solid  or  gas  is  used  to  couple  the 
differential  equations  and  provide  a means  of  solution.  The  equation 
of  state  expressed  in  its  most  general  form  is  {Ref  8:294): 

u2 

1 - f(P,v)  - E - E (4) 

where  I is  internal  energy  and 


Assumptions  and  Resulting  Equations 

Simplifying  assumptions  are  necessary  to  reduce  the  multi-dimen- 
sional equations  describing  shock  propagation  to  the  less  general  one- 
dimensional equations  solved  in  HUFF . The  assumptions  are  as  follows: 

1.  The  local  fluid  velocity  is  directed  along  and  depends 
on  the  Eulerian  one-dimensional  coordinate. 

2.  Adiabatic  conditions  exist  throughout  the  problem. 

3.  An  artificial  positive  viscosity,  2 , can  be  introduced 

to  spread  out  the  density  and  pressure  discontinuities 

that  develop  at  the  shock  front  due  to  the  adiabatic  assump- 
tion. 

4.  The  problem  is  assumed  to  occur  in  free  space  with  no 
gravity  dependence. 

The  differential  equations  in  one-dimensional  form  for  conservation 
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of  mass,  momentum,  and  energy  are? 


|£  + u|£  + p + (-L.  Z 1)u  - 0 
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where  L is  the  geometry  factor  (L  = 1 for  rectangular,  L = 2 for 
cylindrical,  and  L 3 for  spherical). 


Langrangian  Coordinate  System 

The  spatial  mesh  used  in  HUPP  is  constrained  to  conserve  mass  be- 
tween the  mesh  points.  The  volume  between  two  mesh  points  is  defined 
to  be  a cell.  Depending  on  the  geometry,  the  mesh  points  define  cell 
boundaries  that  are.  either  infinite  planes,  infinite  concentric  cylin- 
ders, or  concentric  spheres.  Each  cell  contains  a fixed  mass  for  all 
time  - except  during  rezoning  in  the  blast  problem  where  cells  are  re- 
defined by  removing  every  other  intermediate  cell  boundary.  Requiring 
a cell  to  contain  a fixed  mass  defines  vhat  is  called  a Lagrangian 
mesh.  The  position  of  the  cell  boundaries  are  followed  in  time  and  no 
fluid  moves  through  the  ceil  boundaries  (Ref  5t62). 

In  a Lagrangian  coordinate  system  the  mass  of  each  cell  remains 
fixed.  As  time  progresses,  the  cell  boundary  position  changes.  The 
density  of  a cell  at  time,  t , is  given  by 
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p(x) 
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■ po(xo> 


3x  t 


where  xQ  and  pQ  are  the  original  position  and  density  respectively 
The  velocity  of  a cell  boundary  is  given  by 


(10) 


The  Lagrangian  equations  for  one-dimensional  shock  propagation 
expressing  conservation  of  mass,  momentum,  and  energy  are 


3u  . (L  - l)Pu 

ir— 
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(ID 


3u  3p 
Po3t  + 'Sx^- 


L-l 
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(12) 
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Lagrangian  Mesh 

The  Lagrangian  mesh  material  and  cell  boundary  properties  are  de- 
picted on  Figure  1.  The  cell  boundary  positions  and  velocities  begin 
with  x^  and  at  the  left  wall  of  the  mesh.  The  material  prop- 
erties are  cell  centered  values  with  ' an,a  located 
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Figure  1.  Labeling  Scheme  for  Lagrangian  Cells 
and  Materials  Properties 

between  the  cell  boundaries  at  Xf  and  xi+i.  Program  HUFF  is  dimen- 
sioned to  contain  a grid  of  101  values  (100  cells) . can  be  the  left 

wall  or  center  of  a rectangular  geometry  problem  while  it  usually  rep- 
resents the  center  of  a spherical  or  cylindrical  geometry. 

Finite  Difference  Scheme 

The  choice  of  a properly  centered  finite  difference  scheme  is  crit- 
ical in  developing  a stable  code  to  propagate  the  shock.  The  finite 
difference  scheme  must  approximate  tha  solution  in  a way  that  strictly 
conserves  mass,  momentum,  and  total  energy.  Any  difference  scheme 
that  is  properly  centered  in  both  time  and  space  can  be  expressed  in  a 
conservative  form  (Ref  9:283).  Inadvertent  changes  in  the  material 
conditions  of  any  cell  in  the  mesh  can  lead  to  propagation  of  unwanted 
extraneous  waves  through  the  mesh. 


The  difference  equations  used  in  MUFF  are  those  developed  by 
Richtmyer  and  Morton.  The  equations  are 


n+1  n 

Xa  - Xj 
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P Ax  \(x0)j' 


v n+1  _ 1 ~ (x!|+1)L 


Tn+1  _n  . [p?+:L  + P2 
*j  ~ + -J i 
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where  t is  internal  energy,  Q is  artificial  viscosity,  and  a is 
the  approximate  number  of  cells  over  which  the  shock  discontinuity  is 
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spread  (Ref  8:318) . 

Internal  energy  is  eliminated  from  the  equations  by  coupling  Eq 
(17)  and  Eq  (18) . Equation  (18)  is  the  equation  of  state  for  the  mat- 
erial in  the  problem.  The  ideal  gas  equation  of  state  is  used  for  air 
and  the  Gruneisen  equation  of  state  for  solids . 

Gas . The  finite  difference  form  of  the  ideal  gas  equations  of 
state  is  (Ref  5:3) 


n+1  n+1  n+1 

pj  - <Y  " DPj  Ij 


(20) 


Combining  Eq  (17)  and  (20)  to  eliminate  internal  energy,  gives 
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nf  n.  . n+1 

n 

« n+1 

, n+1 

n 

:1s (vj’  - (vj 

- Vj) 
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where  Y is  gamma  and 


(Y  - 1) 


(22) 


A coinplete  derivation  of  Eq  (21)  is  contained  in  Appendix  D. 

Solid.  The  finite  difference  form  of  the  Gruneisen  equation  of 
state  is  (Ref  5:3) 


pn+l  ■ fp  \n+1  + v P31*1  rin+1  - It  )n+1l 
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where  ys  is  the  Gruneisen  constant  and  c is  the  normal  density  speed 
of  sound  in  the  material. 

Combining  Eq  (17)  and  (23)  to  eliminate  internal  energy  gives 


A complete  derivation  of  Eq  (27)  is  contained  in  Appendix  D. 

The  difference  equations  presented  above  are  solved  explicitly  in 
HUFF'b  Subroutine  HYDRO.  The  solution  is  examined  in  more  detail  in 
Appendix  A. 

Stability  Conditions 

Constraints  must  be  placed  on  the  size  of  the  time  step,  At  , to 


r 


ensure  sufficient  continuity  is  maintained  in  the  solution  of  material 
properties  in  the  mesh,  if  the  time  step  is  too  large,  local  instabil- 
ities will  develop  overwhelming  the  solution. 

The  Lagrangian  stability  condition  provided  by  Richtmyer  and 
Morton  (Ref  8i298) .is  used  to  specify  the  time  step  in  HUFF.  The  equa- 
tion is 


C.At 
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,n 
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<j+l  ‘ xj 
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where  CA  is  the  local  adiabatic  sound  speed  (Ref  5:6),  The  expression 
for  a gas  is 


(29) 


and  for  a solid  is 
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where 
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f A complete  derivation  of  Eq  (29)  and  (30)  is  contained  in  Appendix  D. 

An  additional  stability  condition  is  included  in  HUPP  to  ensure  the 
cell  boundary  positions  do  not  overlap.  This  condition  is 

(Uj+l  " Uj}  At  < ,5(Xj+l  " Xj>  (32) 


This  condition  becomes  & contributing  factor  when  converging  cell  boundary 
velocities  become  very  high  and  density  valuds  in  the  cell  have  not  had 
time  to  respond. 

The  time  step  is  calculated  at  the  beginning  of  subroutine  HYDRO 
fo.v  all  the  cells  in  the  mesh.  The  smallest  time  step  value  is  used  for 
all  mesh  point  calculations  during  the  same  iteration.  The  next  iter- 
ation begins  by  computing  a new  time  step. 


Variable  Gamma  Approximation 

The  actual  value  for  gamma  in  the  gas  equation  of  state  is  not  a 
constant  for  real  gases.  The  value  ranges  from  1.4  under  normal  atmo- 
spheric conditions  in  air  to  a value  as  low  as  1.1  under  high  internal 
energy  conditions.  A variable  gamma  is  used  in  HUFF  to  calculate  the 
time  step  and  the  new  pressures  in  a gas. 

A "semi-physical  fit"  for  gamma  as  a function  of  internal  energy 
for  air  was  developed  by  h.  R.  Doan  and  G.  H.  Nickel  (Ref  2i5) . A 
rough  fit  to  this  function  was  developed  for  use  in  HUFF  calculations. 
The 'HUFF  rough  fit  gamma  function  with  x expressed  in  the  units  1010 
ergs/gm  is 

for  i A .1315 


(Y  - 1)  - .3981 


(33) 


for  .1315  - 1 - 1.0 


(Y  - 1)  - -.0399(1  - .1315) ^ + .3981  (34) 


for  I * i.o 


(Y  - D - .16  + 


.624 
(2  + 1) 


(35) 


A density  correction  is  included  for  values  of  I - 1.0  and  is  given  by 


(Y  - 1)  - (Y  - D(e/p0  )a(I) 


(36) 


where 

a(l)  ■ .0577{IjOG(I)  } - .0218{LOG(I)  }2  + .0035{lOG (I) }3  + 

.0002{bOG(l)}4  (37) 

Figure  2 compares  the  HUFF  fit  to  the  Doan  and  Nickel  equation  of  state 
for  gamma  in  air. 
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Figure  2.  Comparison  of  HUFF  Rough  Fit  to  the  Doan  and  Nickel 
Equation  of  State  for  Gamma  in  Air 


III.  Problem  Definition  and  Solution 


Program  HUFF  was  developed  to  be  uned  with  a wide  variety  of  one- 
dimensional problems.  Each  task  accomplished  by  the  computer  code  is 
contained  in  modular  units  called  subroutines . All  real  or  integer 
values  commo'  mong  subroutines  are  passed  through  the  subroutine  argu- 
ment list.  This  type  of  programming  concept  is  called  modular  pro- 
gramming (Ref  11) . 

Program  HUFF  is  designed  to  solve  any  one-dimenaional  single  ma- 
terial problem  involving  hydrodynamic  shock  propagation  in  a solid  or 
in  air.  A spherical  blast  problem  in  air  and  a solid  compression  problem 
are  included  with  the  program  to  demonstrate  huff's  capabilities.  The 
spherical  blast  problem  will  be  discussed  first  followed  by  the  solid 
compression  problem.  A discussion  of  the  flexibility  of  using  HUFF  in 
solving  additional  hydrodynamic  problems  is  presented  last. 
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Spherical  Blast  in  Air 

The  spherical  blast  problem  is  defined  by  subroutine  BLAST.  The 
initial  conditions  are  based  on  the  followings 

1.  A large  fraction  of  the  yield  is  absorbed  by  the  air  as 
internal  energy, 

2.  All  the  absorbed  energy  is  deposited  before  the  air  has 
a chance  to  respond. 

3.  The  time  required  to  deposite  the  blast  energy  is  defined 
as  the  time  to  hydroseparation  (Ref  4 >65). 

4.  The  absorbed  blast  energy  is  deposited  as  internal  enargy 
in  an  isothermal  sphere  with  a radius  equal  to  the  fireball 
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Table  I 


Average  Atmoshperio  Parameters  for 
Mid-Latitudes  (Ref  4:104) 


Altitude  Scaling  Factors  Altitude  Sealing  Factors 


(feet) 

SP 

SD 

ST 

(feet) 

SP 

SD 

ST 

00 

1.00 

1.00 

1.00 

15,000 

0.56 

1.21 

1.28 

1,000 

0.96 

1.01 

1.02 

20,000 

0.46 

1.30 

1.39 

2,000 

0.93 

1.03 

1.03 

25,000 

0.37 

1.39 

1.53 

3,000 

0.90 

1.04 

1.05 

30,000 

0.30 

1.50 

1.68 

4,000 

0.86 

1,05 

1.07 

35,000 

0.24 

1.62 

1.86 

5,000 

0.83 

1.06 

1.08 

40,000 

0.19 

1.75 

2.02 

10,000 

0,69 

1.13 

1.17 

radius  at  hydroseparation  (Ref  4:65) , 

Subroutine  BLAST  contains  tho  proper  input  values  for  a one  kiloton 
sea  level  explosion.  The  subroutine  scales  the  programmed  values  cal- 
culating the  proper  initial  conditions  for  a variety  of  blast  problems. 

The  yield  of  the  explosive  and  the  altitude  parameters  provided  in  Table 
1 are  input  by  the  user  to  define  the  type  of  blast  problem.  The  param- 
eters are  defined  as 

SP  - (38) 

o 

SD  - (I®) 1/3  (39) 

P 
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(40) 


ST  - (£0)1/3(10,1/2 

p T 

where  pQ  and  TQ  are  normal  sea  level  atmosphere  pressure  and  tem- 
perature and  p and  T are  the  pressure  and  temperature  at  altitude. 
These  parameters  (SP,  SD,  and  ST)  are  used  only  to  define  the  proper 
initial  conditions  for  a variety  of  altitude  and  yield  problems  and 
are  not  used  after  the  hydrodynamic  calculations  begin. 

Compression  of  a Solid 

The  compression  problem  is  defined  by  subroutine  SQUEEZE.  An 
initial  constant  pressure  is  applied  in  the  outside  cell  of  the  mesh. 

The  compression  problem  will  compress  the  solid  in  any  of  the 
three  problem  geometries  (rectangular,  cylindrical,  spherical).  The 
example  problem  presented  in  Chapter  IV  uses  rectangular  geometry. 

The  Gruneisen  equation  of  state  used  in  program  HUFF  limits  the 
amount  of  overpressure  that  can  be  applied  to  the  solid.  The  maximum 
density  that  can  be  achieved  is  less  than 


£ - Ya  + 

po  Y„  - 1 


(41) 


At  this  density,  subroutine  HYDRO  will  calculate  an  infinite  value  for 
pressure . 

• Defining  Other  Problems 

Many  additional  problems  can  be  defined  and  solved  by  program  HUFF 
by  defining  the  initial  and  boundary  conditions  in  a problem  defining 
subroutine  similar  to  BLAST  or  SQUEEZE.  A general  understanding  of  a 
few  basic  concepts  used  in  program  HUFF  will  ba  helpful  whan  defining 
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1, 


a new  problem. 

The  HUFF  program  flow  pattern  1b  depicted  in  Figure  3,.  The  program 
begins  by  calling  subroutine  FLAGS  to  establish  the  "sign  posts"  that 
identify  the  route  the  program  will  follow  in  solving  the  problem;  The 
main  program  then  transfers  complete  control  of  the  problem  solution  to 
subroutine  OONTROL. 

Subroutine  CONTROL  directs  the  routing  of  the  program  during  the 
solution  process.  The  "problem  defining  subroutine"  is  called  before 
beginning  the  solution  to  establish  the  initial  conditions  of  the  problem 
All  the  remaining  task  oriented  subroutines  are  called  at  the  proper 
times  during  the  solution  as  the  program  cycles  through  each  time  step. 

Subroutine  Hydro  calculates  all  the  new  material  properties  of  the 
mesh  for  one  time  step  each  time  it  is  called  by  subroutine  CONTROL. 

The  only  exception  is  the  two  exterior  boundary  velocities  (u-^  ana  u10i)  . 
If  the  boundary  velocities  are  constant  throughout  the  problem  solution, 
no  additional  adjustments  to  the  HYDRO  calculations  are  required.  If 
the  boundary  velocities  change  during  the  problem,  an  additional  entry 
from  subroutine  HYDRO  to  the  "problem  defining  subroutine"  will  be 
necessary  to  calculate  the  outer  boundary  velocities.  This  call  must 
be  accomplished  after  subroutine  HYDRO  calculates  the  time  step  and 
before  the  new  cell  positions  are  confuted. 

The  information  provided  above  should  be  sufficient  to  allow  a pro- 
grammer to  use  HUFF  to  solve  a variety  of  additional  problems  beyond 
those  provided  in  the  original  program.  All  of  the  task  oriented  sub- 
routines contained  in  HUFF  are  well-documented  to  aid  the  user  if 
further  modifications  are  required. 
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Figure  3.  Program  HUFF  Flow  Diagram 


IV.  Results  and  Discussion 


This  chapter,  along  with  Appendices  B and  C,  presents  and  discusses 
the  results  for  two  sample  problems.  The  first  problem  discussed  is  a 
one-kiloton  nuclear  burst  in  infinite  sea  level  air.  The  second  problem 
is  the  one-dimensional  compression  of  uranium  both  in  rectangular  and 
spherical  geometries.  The  results  of  these  problems  demonstrate  the 
usefulness  and  limitations  of  HUFF  as  an  effective  first  order  hydro- 
dynamic  computer  code. 

1_  KT  Explosion 

The  one-kiloton  nuclear  explosion  in  infinite  sea  level  air  was 
chosen  as  a sample  problem  because  of  the  availability  of  other  results 
to  check  those  of  HUFF.  The  time  used  in  the  initial  conditions  was 
.0368  milliseconds  corresponding  to  the  time  of  hydroseparation  (Ref  4: 
65 , 103) . The  radius  of  the  isothermal  sphere  at  hydroseparation  is 
450  cm.  An  input  overpressure  of  1.5  x 1010  dynes/cm2  was  used  to  rep- 
resent the  total  energy  deposited  in  the  isothermal  bubble  by  the  time 
of  hydroseparation.  The  input  pressure  corresponds  to  the  total  energy 
deposited  in  the  air  by  the  time  of  hydroseparation  (Ref  4:64)  . The 
time  calculated  by  HUFF  is  the  time  after  the  nuclear  blast  detonation. 

The  results  from  HUFF  were  compared  with  those  from  the  Nuclear 
Blast  standard  (1KT)  (Ref  7)  and  BEERAY  (Ref  1:36).  The  results  com- 
pared quite  well  with  the  standard  showing  a significant  improvement 
over  the  results  of  BEERAY, 

HUFF  was  compared  with  the  1KT  standard  for  all  the  times  pub- 
lished in  the  standard  over  the  time  range  of  three  millisecond  to  two 
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seconds  after  detonation.  HUFF'S  accuracy  deteriorated  after  two 
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seconds  due  to  its  failure  in  simulating  weak  shock  propagation  at 
great  distances . The  range  during  the  first  two  seconds  is  from  27 
meters  to  2.6  kilometers. 

The  shock  position  calculated  by  HUFF  is  within  three  meters  of 
the  standard  throughout  the  range  from  27  meters  to  beyond  one  kilo- 
meter. At  distances  beyond  one  kilometer  to  over  three  kilometers, 
it  agrees  with  the  standard  within  two  percent.  Almost  all  of  this  un- 
certainty can  be  attributed  to  the  relatively  large  cells  of  the  mesh 
(.5  meters  per  cell  at  27  meters  and  8 meters  per  cell  at  one  kilometer). 

The  maximum  shock  pressure  agrees  generally  within  seven  percent 
of  the  standard  with  a maximum  error  of  13  percent  occurring  around  50 
milliseconds.  The  pressure  is  generally  higher  at  early  times  than  the 
standard  and  lower  than  the  standard  at  late  times.  HUFF • s deviation 
from  the  standard  could  be  the  result  of  neglecting  to  compensate  for 
the  energy  lost  and  gained  in  the  air  by  radiation  from  the  fireball 
and  the  early  shock.  A significant  percentage  of  the  total  weapon 
yield  is  carried  away  by  radiation  and  deposited  ahead  of  the  chock. 

This  energy  is  absorbed  by  the  air  changing  the  normal  conditions 
asflumed  by  HUFF,  Both  pressure  and  density  are  sensitive  to  this 
additional  energy  transfer. 

The  maximum  Bhock  density  agrees  generally  within  five  percent  of 
the  standard.  The  largest  deviation  was  12  percent  occurring  around  50 
milliseconds  corresponding  to  the  high  pressure  present  at  this  time. 
Much  of  this  additional  energy  represented  by  a high  pressure  and  den- 
sity would  be  converted  to  radiation  energy  and  transmitted  ahead  of  the 
shock  due  to  the  high  taqpsratures  present.  HUFF  does  not  contain  a 
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thermal  radiation  code  to  handle  this  condition  properly.  At  later 
times,  pressure  and  dansity  both  return  to  closer  agreement  with  the 
standard  and  fall  below  the  standard  at  great  distances.  This  agrees 
with  what  would  be  expected  from  a code  neglecting  radiation  transport. 

The  maximum  material  velocity  at  the  shock  agrees  within  nix  per- 
cent of  the  standard  throughout  the  strong  shock  propagation  phase. 

\ * 

This  covers  a range  from  27  meters  to  over  two  kilometers.  The  material 
velocity  remains  slightly  above  the  standard  throughout  the  problem. 

In  general,  HUFF  agrees  with  the  Nuclear  Blast  Standard  peak  values 
within  10  percent  with  few  exceptions.  This  is  a remarkable  correlation 
considering  the  small  spatial  mesh  used  (nine  cells  initially  and  50- 
100  cells  at  later  times)  and  the  absence  of  a subroutine  to  handle 
thermal  radiation.  Tables  II  and  III  give  a comparison  of  HJFF's  re- 
sults for  three  milliseconds  and  one  second  respectively . 

HUFF  has  also  maintained  the  efficiency  of  BEERAY  even  though  addi- 
tional steps  were  added  to  improve  HUFF'S  accuracy.  The  addition  of  a 
variable  gamma  in  the  equation  of  state  for  air  significantly  Increased 
the  computational  time  required  to  compute  each  time  step  (cycle) . 

This  additional  computational  time  per  cycle  was  offset  by  the  incorpor- 
ation of  a "hydromatic"  time  step  control  which  reduced  the  number  of 
cycles  required  to  reach  two  seconds  from  9000  iterations  in  BEERAY  to 
2620  in  HU5F.  Figure  4 through  6 further  display  HUFF'S  performance. 


Compression  of  Uranium 

The  rectangular  geometry  compression  problem  was  selected  to  test 
the  equation  of  state  and  the  solution  process  used  in  HUFF  in  solving 
problems  involving  solids.  The  material  selected  for  the  problem 
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Table  XI 


Comparison  of  the  Results  of  HUFF  with 
the  One  Kiloton  Standard 


Time  « 3.00x10" 

’3  sec  (471  Time 

Iterations  for  HUFF) 

1KT  standard 

HUFF 

%Difference 

Shock  Position  (cm) 

2.9xl03 

2. 7x10 3 

-7 

Shock  Overpressure 
(dynes/cm2) 

1.45x10s 

1.52xl08 

6 

Shock  Overdensity 
(gm/cm3) 

7.50xl0“3 

7.44xl0"3 

-1 

Shock  Material  Velocity 
(cm/sec) 

3.45xl05 

3.38xl05 

-2 

(Ref  7:22-26) 

Table  XXI 

Comparison  of  the  Results  of  HUFF  with 
' the  One  Kiloton  Standard 

Time  - 1.000 

sec  (2340  Time 

Iterations  for  HUFF) 

1KT  standard 

HUFF- 

%Difference 

Shock  Position  (cm) 

5.04xl04 

5.05xl04 

0 

Shook  Overpressure 
(dynes/cm2) 

2.3xl04 

2.2x10s 

-4 

Shock  Overdensity 
(gm/cm3) 

2.0xl0"4 

1.94xl0"4 

-3 

Shook  Material  Velocity 
(em/sec) 

5.09xl03 

5. 04x10 3 

-1 

(l»f  7 1 74-81) 
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Distance  (meters) 


Figure  4.  Comparison  of  HUFF's  Ground  Range  vs  Time 
with  the  1KT  Standard 
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Figure  5.  Comparison  of  HUFF'S  Overpressure  vs  Radius 
with  the  1KT  Standard 
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Iterations 


Figure  6.  Comparison  of  HUFF'S  Time  vs.  Iterations 
with  BEERAY  in  the  Blast  Problem  (total  ex- 
ecution time  for  both  HUFF  and  BEERAY 
are  about  equal  due  to  the  additional  time 
per  iteration  required  by  HUFF  to  include 
the  variable  gamma  computations) 
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was  uranium  with  three  percent  molyldenum  by  weight  for  which  well  doc- 
umented data  is  available  on  the  Rankine -Hugoniot  parameters  (Ref  6: 

534) . This  same  material  was  used  in  BQERAY  providing  an  additional 
comparison  for  HUFF'S  results. 

Appendix  C contains  the  input  data  and  some  of  HUFF'S  plotted  re- 
sults for  this  problem.  Figure  7 shows  the  peak  shock  pressure  vs, 
time  and  Figure  8 shows  the  peak  density  vs.  time,  a constant  pressure 
of  one  megabar  was  applied  to  the  outer  cell  of  the  mesh  propaqating  a 
shock  inward.  Figure  9 shows  a plot  of  position  vs.  time  for  the  shock. 

Table  IV  contains  the  results  of  HUFF  compared  to  the  Rankine- 

0 * 

Hugoniot  parameters  and  BEERAY's  results.  The  values  given  for  HUFF 

g 

are  an  average  of  15  data  points  taken  from  a time  range  from  2 x 10 
-5 

to  2 x 10  seconds.  This  corresponds  to  a shock  position  range  from 
9 to  1.5  cm  respectively  in  the  10  cm  initial  radius  sphere  of  uranium. 

The  shock  pressure  at  the  Bhock  front  is  slightly  higher  than  the 
pressure  behind  the  shock.  This  results  in  a slightly  higher  density 
and  material  velocity  at  the  shock  front.  A two  percent  increase  in 
pressure  at  this  point  on  the  Hugoniot  curves  corresponds  to  a density 
increase  of  about  one  percent  due  to  the  slope  of  the  curve  at  this 
point  (Ref  2:87)  . 

Appendix  H contains  the  plotted  results  for  the  HUFF  compression 
problem  in  a spherical  geometry.  In  this  geometry,  the  parameters 
change  as  the  shock  travels  toward  the  center  of  the  sphere.  Comparing 
the  density  and  velocities  at  the  point  in  the  sphere  where  the  pressure 
is  one  megabar  gives  similar  results  to  that  observed  in  the  rectangular 
geometry , 
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Table  IV 


Comparison  of  HUFF'S  Compression  of  Uranium 
with  the  Rankine-Hugoniot  Parameters 
and  BEERAY ' s Results 


Time  - 

1x10“^  sec 

(425  Time 

Iterations  for  HUFF) 

BEERAY 

Standard  HUFF 

%Deviation 

Shock  Pressure 
(megabars) 

1*0 .06 

1 

1.02*0.01 

2 

Shock  Velocity 
(105cm/sec) 

-4.42 

-4.44 

-4.44*0.37 

0 

Shock  Density 
(gm/cm^) 

25.5*0.3 

25.465 

25.73-0 .03 

1 

Particle  Velocity 
(10^cm/sec) 

-1.24*0.021  -1.22 

-1.26*0.01 

3 

(Ref  1 ; 39) 


30 


V.  Conclusions  and  Recommendations 

HUFF  is  a simple,  stable  hydrodynamic  computer  code  based  on  con- 
servation principles.  The  computational  time  is  relatively  short  pro- 
viding results  that  are  reasonably  accurate.  The  conclusions  drawn 
from  the  results  presented  in  Chapter  IV  will  be  listed  first  followed 
by  the  recommendations. 

Conclusions 

The  conclusions  based  on  the  results  of  the  two  sample  problems  are ; 

1.  HUFF  is  a relatively  fast  computer  code  because  the  spatial 
mesh  is  limited  to  100  cells,  the  equations  of  state  are 
simple,  and  a variable  time  step  control  is  employed  to 
optimize  the  time  step.  In  the  blast  problem  rezoning  the 
mesh  at  later  times  results  in  a significant  reduction  in 
the  computation  time  when  propagating  the  shock  to  great 
distances . 

2.  HUFF  provides  reasonably  accurate  results  when  propagating 
strong  shocks . The  results  are  within  13  percent  for  pressure 
and  density,  5 percent  for  material  velocity,  and  2 percent 
for  shock  position  in  the  blast  problem  and  within  2 per- 
cent for  all  parameters  (position,  density,  velocities)  in 
the  compression  problem. 

3.  HUFF  is  a simple  code  due  to  the  techniques  of  nodular  pro- 
gramming used  in  its  development  and  the  reduction  of  input 
data  required  to  operate  the  sample  programs, 

4.  HUFF  is  a flexible  code  capable  of  solving  strong  shock 
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propagation  within  larger  codes  or  solving  additional  problems 
with  minimum  alteration  of  the  existing  code.  This  capability 
is  a result  of  the  task  oriented  independent  subroutines  used 
in  HUFF. 

Recommendations 

The  following  recommendations  present  the  author's  view  concerning 
possible  extensions  to  IttJFF  that  may  improve  its  accuracy  and  provide 
an  even  better  code. 

1.  Develop  additional  subroutines  to  handle  the  radiation  trans- 
port problem  for  early  times  after  a nuclear  detonation  and 
couple  these  subroutines  with  the  blast  expansion  problem. 

2.  Develop  additional  subroutines  to  handle  the  neutronics  in- 
volved in  exploding  the  compressed  uranium  and  couple  these 
subroutines  with  the  compression  problem  in  HUFF,  An  addi- 
tional equation  of  state  may  be  necessary  to  handle  the  ex- 
tremely high  resulting  temperatures  and'  pressures . 

3.  Combine  the  first  two  recommendations  by  coupling  the  energy 
output  from  the  explosion  to  the  expansion  problem  providing 
a complete  computer  code  that  calculates  the  entire  explosion 
sequence . 
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Appendix  A 


The  HUFF  Code 

This  appendix  describes  the  HUFF  computer  code  written  in  the  CDC 
FORTRAN  Extended  4.4  language.  The  first  section  discusses  in  general 
the  program  flow  and  defines  the  specific  task  of  each  subroutine.  The 
final  section  defines  the  primary  program  variables . a complete  listing 
of  program  HUFF  is  contained  in  Appendix  E« 

Program  Flow 

The  HUFF  code  contains  seven  subroutines  and  two  function  subprograms . 
Each  subroutine  performs  a specific  task  in  solving  the  assigned  problem. 
All  variables  common  among  subroutines  are  passed  through  the  appropriate 
argument  lists . 

Main  Program . The  main  program  is  very  short,  its  basic  purpose 
is  to  identify  the  input,  output,  and  local  files  used  by  HUFF  throughout 
the  program.  The  main  program  calls  subroutine  FLAGS  and  then  transfers 
control  cf  the  program  to  subroutine  CONTROL. 

Subroutine  FLAGS.  Subroutine  FLAGS  reads  two  input  data  records 
that  define  the  physical  operation  of  the  program.  These  "flags"  are 
used  during  the  program  to  identify  the  program  flow  sequence. 

Subroutine  CONTROL . Subroutine  CONTROL  controls  the  program  execu- 
tion. All  the  primary  task  oriented  programs  are  called  in  the  proper 
sequence  by  subroutine  CONTROL. 

Subroutine  BLAST.  Subroutine  BLAST  is  the  "problem  defining  sub- 
routine" for  the  blast  problem.  This  subroutine  has  two  entry  points. 

A call  to  blast  will  sat  the  initial  conditions  of  ths  problem.  A call 


to  EXPAND  will  enter  subroutine  BLAST  and  calculate  additional  param- 
eters for  output  from  the  program. 

Subroutine  SQUEEZE . Subroutine  SQUEEZE  is  the  "problem  defining 
subroutine"  for  the  compression  problem.  This  subroutine  also  has  two 
entry  points.  A call  to  SQUEEZE  sets  the  initial  conditions  while  a 
call  to  CRUSH  enters  later  in  the  subroutine  and  calculates  the  same 
parameters  computed  by  EXPAND  in  the  blast  problem. 

Subroutine  HYDRO.  Subroutine  HYDRO  is  the  heart  of  HUFF.,  All  the 
hydrodynamic  shock  propagation  calculations  are  accomplished  by  thiB 
subroutine.  The  subroutine  begins  by  calculating  a proper  time  step 
based  on  the  aurrent  material  properties . HYDRO  then  Calculates  the  new 
material  properties  advancing  the  mesh  through  one  time  step.  The  sub- 
routine contains  an  additional  time  step  control  to  optimize  the  com- 
puted time  step.  This  additional  "hydromatic"  control  varies  the  com** 
puted  time  step  forcing  the  primary  shock  to  cross  a cell  boundary  each 
5-10  cycles  through  subroutine  HYDRO, 

Subroutine  REZONE.  Subroutine  REZONE  extends  the  distance  the 
expanding  blast  shock  can  be  followed  by  rezoning  the  nosh  each  time  the 
shock  reaches  cell  <'4.  This  subroutine  carefully  conservep  the  material 
instabilities  present  before  rezone  so  they  match  those  present  after  re- 
zone. ThiB  enables  the  rezoned  shock  to  continue  its  propagation  with 
minimum  interruption. 

Function  EIA.  Function  EIA  determines  the  proper  internal  energy 
corresponding  to  a given  pressure  and  density.  This  subroutine  is  nec- 
essary since  gamma  is  a function  of  internal  energy  and  internal  energy 
ie  not  directly  computed  and  maintained  by  subroutine  HYDRO.  Subroutine 
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HYDRO  and  subroutine  REZONE  both  use  function  El A. 

Function  VGAMMA.  Function  VGAMMA  contains  the  HUFF  rough  fit  for 
the  variable  gamma  function  used  in  the  equation  of  state  for  air. 

This  function  determines  gamma  based  on  a given  internal  energy  and 
density.  Function  VGAMMA  is  used  in  both  subroutines  HYDRO  and  REZONE . 

Subroutine  COPY.  Subroutine  COPY  handles  the  recording  of  results 
for  the  program.  A complete  listing  of  all  the  program  arrays  are 
printed  on  a local  file  called  files  for  each  result  time  requested  by 
the  user.  A rough  plot  of  the  requested  results  is  also  made  with 
appropriate  calls  to  PLOTSM. 

Subroutine  PLOTSM.  The  last  subroutine  contained  in  HUFF  is  sub- 
routine PLOTSM.  This  subroutine  provides  a rough  plot  of  the  program 
arrays  on  the  output  file . A rough  plot  of  all  the  requested  results 
for  aid  in  analysis  of  the  output  is  provided. 

Obtaining  Formal  Output.  When  the  program  execution  is  complete, 
the  requested  results  are  contained  on  two  local  fileB.  A file  named 
files  contains  the  requested  information.  A file  named  Review  contains 
summary  information.  The  information  can  bs  read  from  these  files  and 
minipulated  into  any  desired  form.  An  additional  program  called  Results 
is  listed  in  Appendix  E that  will  read  the  data  contained  on  the  local 
files  and  request  the  appropriate  formal  plots. 

Main  program  Variables 

This  section  contains  a listing  of  all  the  program  variables  used 
in  HUFF.  Many  of  the  acronyms  suggest  their  specific  use  in  the  program. 
Unless  otherwise  Indicated,  all  program  variables  are  in  the  gram-centi- 
meter-second  system  of  units.  The  program  variables  are i 
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A: 


Temporary  storage  variable. 
Adiabatic  sound  speed  in  a solid. 


C: 

CA: 

CK ( 1) /CYCLE/ICY CLE • 

CK(2)/TIME  s 

CK(3)/DT- 
CK(4)/P0S  s 
CK(5)/OP: 

CK(6)/OD: 

CK(7)/VMAX: 

D(I)  s 
PMAX: 

DTJ; 

DX: 

EI/E1E;  . 

EKJ1/EKJ2/EKI1 : 

EL: 

ETJ1/ETJ2/ETI1 : 
PLAG(2)/IF2j 

■I 

FLAG (3)/IF3: 

FLAG  (4): 

GAMMA: 

GMl/GMONE : 


The  square  of  the  number  of  cells  over 
which  the  shock  discontinuity  is  spread 
(a2  in  Eq  (19)) . 

Counters  for  the  number  of  cycles  through 
subroutine  HYDRO. 

Keeps  track  of  total  real  time  (sec) 
during  the  solution. 

The  time  step  value  (sec)  . 

Shock  position. 

Peak  shock  overpressure. 

Peak  shock  overdensity. 

Peak  material  velocity  at  the  shock. 

Cell  centered  density  array. 

Maximum  cell  centered  density  at  the  shock 

Temporary  local  time  step  storage  variable 

Original  mesh  spacing  of  the  array. 

Internal  energy  in  ergS/gram  and  1010 
ergs/gram  respectively. 

Kenetic  energy  storage  variables. 

Natural  logrithm  of  EIE. 

Total  energy  storage  variables . 

Flag  for  the  number  of  results  times 
requested. 

Flag  for  the  maximum  number  of  allowed 
cycles  through  subroutine  HYDRO. 

Flag  used  to  identify  the  type  of  problem. 

Storage  variable  for  the  Gruneisen  con- 
stant or  gamma. 

Storage  variable  for  gamma  minus  one. 


HIOLD/HINEW: 

HPOLD/HPNEWt 

IALPHA/L i 
ICOUNT/KFIFTY • 
IP/Nl/USAVE  s 
MN/M: 

MPlt 

NM5/NP5 : 

OPOSj 
P(I)  : 

P2/PX  j 
PRESS : 

P'JD  (I , J)  t 

PMAX; 

PNORMs 

PBLAST : 

Q(I)  : 

QMAXi 

RHOZi 
RFO  i 
RADIUS i 

St 


Temporary  storage  variables  used  for  X 
from  the  solid  equation  of  state.  h 

Temporary  storage  variables  used  for  PH 
from  the  solid  equation  of  state. 

Problem  geometry  factors. 

Loop  counters. 

Shook  grid  position. 

Used  as  a variable  length  for  do  loops  in 
subroutine  HYDRO. 

M plus  one. 

HI  minus  5 and  plus  5 respectively. 

Old  position  of  the  shock. 

Material  property  array  for  cell  centered 
pressures . 

Temporary  pressure  storage  variables. 

Input  pressure  for  the  compression  problem. 

Temporary  storage  array  for  cell  centered 
pressure,  velocity,  and  density  respectively. 

Maximum  pressure  at  the  shock  front. 

Normal  atmospheric  pressure  at  the  burst 
altitude . 

Input  pressure  for  the  blast  problem. 

Material  property  array  for  the  artificial 
cell  centered  artificial  viscosity. 

Maximum  viscosity  at  the  shook  front. 

This  occurs  at  the  location  of  the  shock 
by  definition. 

Ratio  of  D(I)/RHO. 

Normal  material  density. 

initial  radius  of  the  solid  in  the  compression 
problem. 

Temporary  storage  variables  for  various 
functions  of  gamma. 
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SP/SD/ST: 

Altitude  scaling  parameters  for  pressure, 
density,  and  time  respectively. 

TSTOPt 

Real  time  limit  for  the  problem. 

TP RAC i 

A variable  fraction  used  to  place  additional 
control  on  the  size  of  the  time  step  to 
increase  program  efficiency. 

TIMITt 

Used  as  a loop  timer  in  calculating  shock 
speeds . 

0(1)  « 

Material  property  array  for  the  cell 
boundary  velocities . 

UJ1/UJ2/UI1: 

Temporary  storage  variables  for  cell 
centered  velocity  in  subroutine  REZONE. 

VG (I)  • 

Variable  gamma  storage  array  in  subroutine 
HYDRO. 

VEL(I) s 

Additional  array  used  to  temporarily  store 
velocity  in  subroutine  HYDRO. 

VRHO» 

The  normal  specific  density  of  the  material. 

VOID/’/NEW  j 

Temporary  storage  variables  for  cell 
centered  specific  densities  used  in  sub- 
routine HYDRO. 

VDOTi 

Specific  density  change  in  one  time  step. 

VSVOLD/VSVNEW  s 

Temporary  storage  variable  used  with  the 
solid  equation  of  state. 

X(I)  : 

Material  property  array  for  the  cell  boundary 
positions . 

XI  t 

Real  variable  equal  to  the  integer  X. 

XMJ1/XMJ2/XMI1 s 

Temporary  storage  arrays  for  the  mass  of 
the  two  old  cells  and  the  new  cell  respec- 
tively used  in  subroutine  REZONE. 

XMOMJ1/XMOMJ2 /XMOMIl t 

Temporary  storage  arrays  for  the  momentum 
of  the  two  old  cells  and  the  new  cell 
respectively  used  in  subroutine  REZONE . 

XNORM: 

A variable  used  in  a scaling  function  rep- 
resenting the  normal  value  of  the  property 
to  be  Scaled. 
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Appendix  B 

Sample  Problem,  1KT  Explosion 

The  data  for  HUFF  is  always  read  by  a list  directed  format  to  sim- 
plify the  input  cards. 

The  input  cards  for  tills  problem  are  as  follows  where  the  items 
in  parentheses  are  not  part  of  the  data  deck: 

10  3000  1 5.0 

.003  .006  ,01  .02  .05  .1  .2  .5  1.  2. 

1.  1.  1.  1. 

Figures  4 through  6 in  Chapter  IV  of  the  text  and  Figures  10  through 
18  of  this  appendix  represent  sample  output  for  this  problem. 


(Card  1) 
(Card  2) 
(Card  3) 
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Figure  10.  Initial  Pressure  Condition  for  1KT  Burst 


0£NSITT(C/CH-»J>  VS  »OSI TIOVt SH>  »T  TIHE  • 1.002E»00  SEC 


Figure  17-  Material  Velocity  for  1KT  Burst  (lsec) 
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Sample  Problem , one  Megabar  Compression  of  Uranium 

The  input  cards  for  this  problem  in  rectangular  geometry  are  as 
follows  where  the  items  in  parentheses  are  not  part  of  the  data  deck: 

xO  4000  2 1.0  (Card  1) 

l.E-5  2.E-5  3.E-5  4.E-5  5.E-5  6.E-5  7.E-5  8.E-5 

9.E-5  l.E-4  (Card  2) 

1 2.03  18.45  2.56E5  10.  1.E12  (Card  3) 

The  data  1b  in  list  directed  format  to  simplify  defining  the  problem. 

For  a spherical  compression  problem  change  Card  3 to  the  following: 

3 2.02  18.45  2.56E5  10.  1.E10  (Card  3) 

The  compression  sample  problems  maintain  a constant  applied  pressure 
in  the  outer  cell  of  the  mesh  throughout  the  problem  solution.  This  can 
be  changed  to  an  impulse  pressure  by  removing  the  card  that  redefines 
the  outer  cell  pressure  immediately  after  ENTRY  CRUSH  in  subroutine 
SQUEEZE.  A comment  card  is  included  in  the  program  to  identify  the 
location  of  this  card. 

Figures  7 through  9 in  Chapter  IV  of  the  text  aftd  Figures  19  through 
23  of  this  appendix  represent  sample  output  for  the  rectangular  geometry 
problem.  Figures  24  through  32  of  this  appendix  represent  sample  output 
for  the  case  of  spherical  geometry . 
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Figure  23,.  Material  Velocity  Rectangular  Compression  of  Uranium  (lxl0“5sec) 


Figure  22.  Artificial  Viscosity  for  Rectangular 
Comp res sion  of  Uranium  (lxl0~5sec) 
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Figure  23.  Time  vs  Iterations  for  Rectangular 
Compression  of  Uranium 
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Figure  27.  Density  for  Spherical  Compression  of  Uranium 


0 100  200  300  400  500  600  700  800  900 


Ite rations 


Figure  32.  Time  vs  Iterations  for  Spherical 
Compression  of  Uranium 
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Derivations 

The  derivations  presented  in  this  appendix  develop  the  key  equa- 
tions used  in  program  HUFF.  Coupling  of  the  equation  of  state  with  the 
Richtmyer  and  Morton  difference  equation  for  energy  will  be  presented 
first.  Next,  the  derivation  of  the  adiabatic  sound  speed  is  presented. 
The  last  derivation  describes  the  basis  for  the  equations  used  in  re- 
zoning  the  mesh  in  the  blast  expansion  problem. 

Equation  of  state  Coupling 

Hydrodynamic  shock  propagation  results  from  on  imbalance  of  the 
interrelated  material  properties  of  a gas  or  solid.  Each  particle  inter- 
acts with  its  adjacent  particle  in  an  effort  to  maintain  or  reestablish 
an  equalibrium  condition  in  the  material.  This  interaction  between 
particles  is  simulated  in  Huff  on  a macroscopic  level  by  packaging  units 
of  material  in  Iiagrangian  cells.  Each  cell  interacts  with  its  adjacent 
cell  to  maintain  or  establish  equalibrium  throughout  the  mesh. 

The  solution  of  the  conservation  equations  of  mass,  momentum,  and 
energy  is  dependent  on  the  equation  of  state  of  the  material . The  three 
equations  expressing  the  conservation  laws  contain  four  dependent  vari- 
ables (Ref  Equations  1-3).  The  equation  of  state  defines  an  additional 
relation  between  these  variables  providing  a unique  solution. 

‘ solving  the  conservation  equations  can  be  simplified  by  replacing 
the  expression  for  energy  in  the  energy  conservation  equation  by  its 
equivalence  from  the  equation  of  state  of  the  material.  This  reduces 
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the  problem  to  a solution  of  three  equations  and  three  unknowns  and 
provides  for  a unique  solution. 

The  difference  scheme  used  in  HUFF  provides  an  explicit  solution 
of  the  conservation  equations.  The  energy  conservation  equation  is 
expressed  as 

r 


n+1  n 

Ij  - Ij  + 


Pn+1  + Pn 

pj  + Pj  n+1 

+ Qu 

2 3 


(Vj  - Vj)  - 0 


(17) 


The  next  two  sections  eliminate  internal  energy  using  the  equation  of 
state  for  a gas  and  for  a solid  respectively. 

Gas.  The  equation  of  state  for  a gas  is 


n_n 


(18) 


Substituting  the  egression  for  I in  equation  (20)  into  equation 
(17)  gives 


0n+l 


(Y  - l)ej+1 


(Y  - l)f 


n+1 

n r>  + p11  't 

*-7  t J — 3 + %+i  y 

l)f"  L 2 3^3 


n+1  - Vj) 


« 0 


(42) 


Multiplying  by  2 and  combining  terms  with  p and  Q gives 


l - $ ♦ * 


fir'  ] 


pj  [(vj+1  - vj>  - Vri3  ^ 1 + B3+1(vTl  - vj}  " 0 (43) 
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n+1 

Solving  for  P^  gives 


_n+l 

pj 


n . n 

pi  tsvi 


- <v;+1  - v;»  - 2b5+1<v5+1  - v^, 


sv^1  ♦ <v”+1  - v") 


(21) 


where  S m 


=rnr  • 


Solid . The  equation  of  state  for  a solid  is 


- »B> J + VSf!)  [:j  - 


« )“ 

uh'  j 


(I  )"  . i 

«)  2 


C (v0  - V ) 

{v  - s(v  -V?)}2 

o o 3 


" c(VQ  - Vj>  1 

/o  - S<Vo  ■ *j>Jj 


(23) 


(24) 


(25) 


<YS  + 1) 


(26) 


Solving  for  1^  in  equation  (23)  gives 

n , .n  . Vj  f n . „n  1 

" 'Vj  * f LPj  - PK  j J 


(44) 


Substituting  equation  (44)  into  equation  (17)  gives 
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. Ts 


t - «vr  + <»•;  - v? * «»>;;  j 

* flf  ](vrl-vS"0 


(45) 


Multiplying  by  2ys  and  collecting  terms  with  P and  Q gives 


pn+i^yn+i  + Ys(vn+1  + v",  ] + P;[ys(v5+1  - v“>  - 2V^j-  2(PH)J+1v” 


n+1 


J(Ph)Jv5  V 2Ys[(Ih)"+1  - «H) “ ] + 2Vse^l(v"+1  - V”)  - 0 


(46) 


Solving  for  p"+1  gives 


Pn+1  - fJ<2v"  - YS(VJ+1  - V*j)}  - 2Ys((Ih)"+1  - (1h)J) 

2VJ+1  + YS(V;+1-V") 

+ 2Vj+1XPH)"+1  - 2V^(Ph)"  - 2YSQ5+1(V^+1  - vj) 

2vj+1  + Ys (v j+1  “ Vj)  (27) 


Derivation  of  Adiabatic  Sound  Speed 


The  adiabatic  sound  speed  is  given  by 

V2 


Derivation  of  the  adiabatic  sound  Bpeod  is  presented  first  for  a gas 
and  then  for  a solid. 

Gas. 

P ■ (Y  - l)pl  (20) 


The  derivation  of  pressure  with  respect  to  density  is  given  by 

||-(Y  - 1)1  + (Y  - 1)P§~  (47) 

Replacing  I in  equation  (47)  with  its  equivalence  from  equation 
(20)  gives 


(48) 


Under  adiabatic  conditions  — is  given  by  (Ref  3:5) 

3p 


3l  . P 
* T 
P 


3P 


(49) 


Substituting  equation  (49)  into  equation  (48)  gives 


9P 

ap 


yp 

p 


Thus 


Solid . The  equation  of  state  for  a solid  is 
P - PH  + I*  (I  - IH) 


(50) 


(29) 


(23) 
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(55) 


31, 


3P 


C2V2V02  - c2V3V0 
{VQ  - S(VQ  - V)}3 


For  adiabatic  conditions  — is  given  by 


31 

?P 


PV^ 


(49) 


Substituting  the  above  expressions  and  the  expression  for  I from 
the  equation  of  state  into  equation  (51)  gives 


3P  c2V2Vo(l  + S)  - c2V2S  y£ 
3p  {V0  - s(vo  - V)}3  V 


PV2  “/c2V2V2  - 02V3Vo 


{Vo  - S(VC  - V)  }‘ 


p - p* 


V 


+ - XH 


(56) 


Substituting  the  value  for  PH  into  equation  (56) , collecting 
terms  and  simplifying  gives 


3P  c2V2Vq(1  + S + Ys)  - C2V3S  - Yge2W02 


3P 


{v„  - S(Vn  - V)}J 


c (yo  “ V>V 

{»:~v.  -v)}2  + pv(Tsti’ 


(57) 


Substituting  Yg  " 2S  “ 1 into  equation  (57)  , combining  terms  and 
simplifying  Hives 


3p 

3p 


2Sc2V2V  - Sc2W  2 
2SVP  + ° °_ 

{vo  “ s(vo  “ v>>3 


(58) 
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Thus  the  adiabatic  sound  speed  in  solids  is  given  by 


[2Sc 2V2Vo  - Sc 

2SVP  ♦ 

vo  - s(vo  - V 

p.ezone  Computations  from  Conservation  Requirements 

Rezoning  during  the  blast  expansion  problem  provides  the  capability 
to  use  small  mesh  spacing  during  the  early  phases  of  shock  development 
and  also  follow  the  shock  to  great  distances  later  in  the  problem.  Re- 
zoning is  accomplished  when  the  shock  reaches  the  94th  cell  of  the  mesh. 

Rezonin  compresses  the  information  stored  in  the  100  cell  mesh 
into  the  first  50  *.ellr.t  and  places  undisturbed  air  in  the  last  50  cells. 
The  rezone  process  carefully  conserves  mass,  momentum,  and  total  energy 
to  preserve  the  material  imbalance  present  in  the  mesh  prior  to  rezoning. 

This  section  presents  the  derivations  that  provide  the  equations 
used  in  rezoning.  The  derivation  begins  with  defining  the  location  of 
the  new  cells  and  proceeds  in  the  order  of  solution  used  in  subroutine 
REZONE. 

Position . The  new  cell  locations  in  the  mesh  are  formed  by  removing 
the  even  numbered  cell  boundaries  and  combining  the  adjacent  cells.  The 
old  cell  positions  use  the  subscript  j and  the  new  cell  positions  use 
tiie  subscript  i . The  location  of  each  new  cell  boundary  is  given  by 

Xi  - X1^ 

(60) 

where  j are  the  even  cell  boundaries  from  2 to  98.  Cell  boundary 
locations  after  cell  48  are  computed  from  the  previous  cell  position 

by 


2W  2 
0 


1/2 


(30) 
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xi  “ Xi-1 


+ DX 


(61) 


where  DX  is  twice  the  original  cell  spacing  existing  before  the  rezone. 

Conservation  of  Mass . Conservation  of  mass  provides  an  explicit 
solution  of  the  new  cell  density  and  mass . The  equation  for  conser- 
vation of  mass  is 


pi<Xi+l  - Xi>  - Vxj+1  " Xj>  + Pj+l(Xj+2  - xj+l> 

The  egression  for  density  is 


Vxj+i " xj' 


- x?+1> 


<*j+2  - *j> 


(62) 


(63) 


where  x^+2  ° x^+1  and  xj  “ xi  • 

Cell  Centered  tfe  loclty . The  time  step  calculations  accomplished 
by  program  HUFF  use  the  cell  boundary  velocities  to  change  the  cell 
positions.  To  conserve  momentum  and  energy  in  the  cells  during  rezone, 
the  average  velocity  of  each  cell  is  required.  The  average  cell  velo- 
city is  computed  by  averaging  the  cell  boundary  velocities  asstiming  the 
velocity  varies  linearly  across  the  cell.  The  average  cell  velocity  is 
given  by 


uj  - uj  + “j+l 
2 


(64) 


* tih 

where  u^  is  the  cell  centered  velocity  for  the  j cell. 

Total  Energy.  The  total  energy  of  each  cell  is  computed  from  the 
cell  centered  velocity  and  the  equation  of  state.  The  equation  for 
total  energy  is  given  by 
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2 


I j 


© 


Si! 


j 


<Y  - l)Pj 


(65) 


where  Ej  is  the  total  energy  of  cell  j and  a variable  gamma  is  used 
confuted  from  the  pressure  of  cell  j . 

Conservation  of  Energy.  Conservation  of  energy  provides  an  expli- 
cit solution  for  the  total  energy  in  the  new  cell.  The  equation  for 
conservation  of  energy  is 


- *Z>Ei 


j? 

'j '^j+1 


p - x-)Ej  + - x:^>e 


3 - X3 

j+l'"j+2  Xj+l#“j+l 


(66) 


The  egression  for  the  new  cell  energy  is 

B<  . pi {xj+i  - xi’Ej  + ej+i«*?*2  - xj+i> 


pi (X j+2 


(67) 


Conservation  of  Momentum.  Conservation  of  momentum  provides  an 
explicit  solution  for  the  cell  centered  velocity  of  the  new  cell.  The 
equation  for  conservation  of  momentum  is 


- 4>~°i  - - *j>«j  + fj+1<xj,2  - 


(68) 


The  expression  for  the  cell  centered  velocity  of  the  new  cell  is 


- - - x3>uj  * - ’’j+l’Vl 

1 »l<*j+2  - xj»  _ 


(69) 


Cell  Boundary  Velocities . lhe  new  cell  boundary  velocities  are 
recomputed  using  a procedure  similar  to  that  used  in  computing  the  cell 
centered  velocities.  The  cell  boundary  velocity  is  computed  by  averaging 
the  cell  centered  velocities  of  the  adjacent  cells  assuming  the  velocity 


i 
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varies  between  the  adjacent  cell  centers  in  a manner  that  gives  an  average 
velocity  at  the  common  boundary.  The  scpression  for  the  cell  boundary 
velocity  of  the  new  cell  is  given  by 


ui-l  + ui 


(70) 


The  cell  boundary  velocity  at  the  center  position,  , is  defined  as 
zero  and  is  not  computed  by  equation  (70) . 

internal  Energy.  The  internal  energy  of  the  new  cell  is  calculated 
by  the  application  of  the  conservation  of  total  energy.  The  kenetic 
energy  of  the  new  cell  is  calculated  first  and  is  given  by 


“i  " 7"  (?i> 

where  KE^  is  the  kenetic  energy  of  the  new  cell  i . Total  energy 
conservation  expressed  in  its  most  simple  form  is 

EjL  - KE^  + 1±  (72) 

where  X^.  is  the  internal  energy  of  the  new  cell. 

The  expression  for  the  internal  energy  of  the  new  cell  is  given 

by, 


- E±  - KEa  (73) 

Pressure . Pressure  is  computed  using  the  equation  of  state.  The 
expression  for  pressure  is 


P 


i " 


(Y  - D&ii 


(20) 


where  (y  - 1)  is  a function  of  internal  energy.  The  value  of  (y  - 1) 
is  computed  directly  using  the  internal  energy  already  computed  for  the 
new  cell. 

Computing  the  pressure  completes  the  rezoning  process  and  the  re- 
zoned  mesh  is  returned  to  subroutine  CONTROL  to  continue  the  Bhock  prop- 
agation. The  mesh  will  again  return  to  subroutine  REZONE  when  the  shock 
position  again  reaches  the  94th  cell. 


Appendix  E 


Program  Listing  of  HUFF 


180s  PROGRAM  HUFF  ? INPUT = /80 . OUTPUT » TAPt 1 = INPUT . TAPE2= OUTPUT 

118=  1,FILES»REVIEN,TAPE7=FILES,TAPF3=REV1EH) 

128=  DIMENSION  FLAC<4I,TABLE(10I 

138=  INTEGER  FLAG 

148=  CALL  FLAG$(FLAG, TABLE, TSTOP) 

158=  CALL  CONTROL. ! FLAG  i T ABLE  i T STOP  I 

148=  PRINT*, "TIKE  LIMIT  EXCEEDED" 

178=  STOP 

188=  END 

190=C 

208=  SUBROUTINE  FLAGS ( FLAG . TABLE , TSTOP ) 

218=  DIMENSION  FLAGI4) , TABLE! 10) 

220=  INTEGER  FLAG 

230-C INITIALIZE  CONTROL  PARAMETERS  

.240=  00  10  1=1,3 

250=  10  FLAG !I 1=0 

260=  DO  30  1=1,10 

270=  38  TABLE (11=0 

288=C  I READ  FLAGS  AND  REAL  TIME  LIMIT  I 

290=C* I FLACi2l=f OF  TABLES  REQUESTED  I 

300=C  I FLAG (31 =f  OF  CTCLES  THROUGH  HTDRO  I 

310=C  I FLAG  141 =TVPE  PROBLEM  il=ELAST,2=SQUEEZEII 

320=  READ*, (FLAG(I) ,1=2,41 , TSTOP 

330*C READ  REQUESTED  TABLE  TIMES 

340=  IF2=FLAC(2) 

350=  READ*, (TABLE(l), 1=1, IFZ) 

360=C UNITE  OUT  INFORMATION  PROVIDED  BT  FLAGS 

370=  IF(FLAC(4).EQ.1)WRITE(2,400| 

360=  1F(FLAG(4I.EQ. 21  WRITE (2, 5H0I 

398=  IF(FLAG(4).NE.1.AND.FLAG(4).NE.2)WRITE(2»600)FLAG(4! 

408=  HRITE(2,200IFLAG(2),(TABLE(I),I=1,IF2! 

410=  URITE(2,300)TSTOP,FLAG(3) 

420=  200  FORMAT!////, T10, "TABLES  FOR  THE  FOLLOWING  ",!3, 

430*  1"  TIKES  HERE  REQUESTED!", //,2<T15,5UPE10.3,2X),/)) 

440=  300  FORMAT (//// ,T10,"hflX  TIKE  10  RUN  IS  ",F4.1,"  SEC." 

. 450=  1,/,T10,"MAX  CTCLES  THROUGH  TIKESTE?  IS", 19,".") 

460=  400  FORHAT(1HI,T40,"THIS  IS  A BLAST  PROGLEM") 

470=  500  FORMAT (1HI.T4D»“THIS  IS  A SQUEEZE  PROBLEM") 

460*  600  FORMAT <1H1,T40,"H0  SUCH  PROBLEM  NUMBER  * ",I3) 

490*  RETURN 


528=  SUBROUTINE  CONTROL  I FLAG i T ABLE i T ST  OP ) 

538=  DIMENSION  X ( 10 1 ! »U ( 10 1 ) t P <13 1 > » D ( 1 01 ) 

548=  1 i TABLE ( 10 ) i CK<7) 

558=  INTEGER  FLAG 

560*C CALL  TO  SET  INITIAL  CONDITIONS  ----- - 

578=  IF(FLAG(4hEQ.l) 

588=  1CALL  BLAST (CA.IALPHAiGAHNAiNl.H.DIiRHO.PtDiIiQ.UiDTiCK) 

570=  IF (FLAG (A) .EG. 21 

600=  1 CALL  SQUEEZE (CAt IALPHAi GAMHAtNl >DX*RHOt iQ>Ui C2iDT»CK) 

m-  N1=0 

620=C LIST  INITIAL  CONDITIONS 

630=  TIHE=CK(2) 

640=  CALL  ICOPT(XiUiPiDf9iFLAGl2bTlHEiTABLEiCKl 

650=  IF3=FLAC(3I 

660=  DO  15  ICYCLE-1 »1F3 

670=C CALL  TO  ADVANCE  ONE  T1HESTEP 

680=  CALL  KYDROlOiPt KiSiUiDXiCAt lALFKAtRKGtCAMllAiNi iDTtTIMEtFLACtCZ) 

690=C CALL  TO  SET  BQUNDART  CONDITIONS  AND  TRACK  SHOCK 

780=  IF(FLAC(4).E9.1i 

710=  1CALL  EXPAMOlCAi IALPHA»GAMHA»N1 iM»DX  »RHOtPtO>XtQfUiDT »CK) 

720=  IF(FLAG(4).E0.2I 

730=  1 CALL  CRUSH4CA? IALPHA»GAHHA»N1  »DKBRH0»P»D*lt»Q»U»C2f DT»CK) 

748=  TIME=CK(2) 

750=C CALL  TO  COPT  REGUESTED  TABLES  

760=  IFITIME.GE.TSTOP.QR.TIHE.GE.TAElEd)) 

770=  1CALL  COP t<XfU»P*0»Qi FLAG ( Z ) »TIMEiTABLEiCK) 

780=C ABORT  IF  TINE  LIMIT  REACHED - 

790=  IF (TIME .CE .TSTOP) RETURN 

800=C CALL  FOR  REZOME  IF  REQUIRED 

810=  IFIN1.LT. 94. OR. FLAG14I  .NE.1ICQ  TO  15 

820=  CALL  REZONE4XiUtP»DiQiDX iGAHMAt ICTCLEtTlHEtNl tRHO) 

830=  15  CONTINUE 

840=C CALL  RESULTS  AND  ABORT  PROGRAM  IF  CTCLE  LIMIT  REACHED  

850=  ' CALL  COPY <X»U*P»D»QiFLAC(2J tTIMEtTABLE»CKl 

860=  PRINTVTQU  REACHED  THE  CTCLE  LIMIT  OF  "iICTCLE 

870=  PRINTLINE  IS  "iTIMEi"  SEC." 

880=  RETURN 

890=  END 

900=r 

910**  SUBROUTINE  HTDROIDiPiIiOiUtDlUCAiLiRHQ»GAHMA»KttDTfTIMErFLAGtC2) 

920=  DIMENSION  D(181)iPU01))I(10U>QU01)iUU01)tFLAC(4) 

930=  ltVC(101)»VEL<ieU 

940=  INTEGER  FLAG 

950=C INITIALIZE  LOCAL  PARAMETERS  

968=  VEL(tl=0 

978=  DT=1. 

eo  a,  MttlUA 

990=  IF(FLAG(4).EO.21H=10I 

1880=  IF(M.LT.15»M>15 


1010*C INITIALIZE  HTDR0WAT1C  TIME  STEP 

1020s  IF1H1.NE.0IGO  TO  10 

1130s  DO  5 1=1,101 

1040=  5 VC(I)=l,4 

1050*  TFRAC8 .1 

1060*  !F<FLAC14).EQ.2iNl*95 

1070*  IF (FLAG { 4) ,E8.1)N!=10 

1060*  ICOUNT*0 

1090*  10  I COUNT *1 COUNT* 1 

1100*  NSAVE'Nl 

tii0*  15  CONTINUE 

1120*  DO  20  0*2 tH 

1130*  lFiFLAC(4).E8.2)G0  TO  17 

1140=0 COMPUTE  CAS  T1MESTEP  PARAMETER 

1150*  A*D(J)/<VG(JUP(J)) 

1160*  GO  TO  18 

1170*  17  CONTINUE 

U80*C COMPUTE  SOLID  T1HFSTEP  PARAMETER  - 

1190*  V*1./D(J» 

1200=  VRH0*1./RHR 

1210*  S*(GAMMA+U72. 

1220*  A*l.  / (2.  *S*P  i J1  »V+  (C2*V«2*VRH0*  ll.+S)  -S* 

1230=  1C2*V*VRH0**Z1 / 1VRH0-S+ IVKHQrV) )**31 

1240*  18  CONTINUE 

1250*0 COMPUTE  LOCAL  T1MESTEP 

1260*  IFIA.LE.0.1CO  TO  20 

1270*  DTJ* (X LM) -I (J) ) «TFRAC*SQRT (A) 

1280*C LOCATE  MINIMUM  TIHESTEF  

1290*  IF(DTJ.LT.DT)DT*DTJ 

1300*  20  CONTINUE 

1313*  DO  25  I*2.M 

1320*C COMPUTE  NEW  CELL  BOUNDARY  VELOCITIES 

1330*  XI*I 

1360*  VELUl*U(mUT»{P(I-l)-P(mQU-l)-Q(I))/ 

1350*  l<RHO«DXH(M)/UXHUDim»<l-H 

1360=C CHECK  CONVERGENCE  VELOCITY  OF  CELL  BOUNDARIES  ----- 

1370*  IF«V£UI-U-VELa!liDT.LT..5MX<I»-X<I-m)G0  TO  25 

1300*  TFRAC*TFRAC*,5 

1394*  PRINT*,"  CONVERGENCE  PROBLEM i TFRAC3** *TFRAC 

1400*  CO  TO  15 

1410*  25  CONTINUE 

1420*  DO  30  I*2»M 

1430*  U(Il*VELtn 

144I*C— COMPUTE  NEW  CELL  BOUNDARY  POSITIONS  

1450*  I(U*XUHDT»UU) 

1461*  30  CONTINUE 
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1470*C COMPUTE  NEW  DENSITIES  AND  SPECIFIC  VOLUMES  

1480'  VRHOM./RHO 

1490s  OMSK =10. 

1500s  DO  40  1-iiH 

1510s  11=1 

1520s  VOIDM./BU) 

1530s  D(I)sRHOi((lliD()««L-((lI-l.)iDnoU/ 

1540s  Kia4iH«L-im«*L) 

1550s  VNEW=1./D(I) 

1560s  VDOT=VNEW-VOLD 

1570SC COMPUTE  NEW  VISCOSITIES  

1580s  Q(ns2.#CA/(VNEW4V0LDU(U(l4l)-U(m«2 

1590=  iF(ua4n-ij(i).GE.0.io(n=0. 

1600=  IF(FIAG(4).EQ.1)G0  TO  36 

1610=C COMPUTE  PRESSURES  FOR  SOLID 

1620=  $s(GAMHA+l.l/2. 

1630=  VSVCLD=VRHO-S4{VRHO-VOLOI 

1640=  VSVNEW=VRHO-Si (VRHO- VMEW1 

1650=  HP0LD=C2*'VRH0-V0L0)/VSV0LDi«2 

1660=  HPNEW=C24  (VRHO- VNEW)  /V$VNEN«2 

1670=  HIOLD=C24((VRHO-VOLD)/VSVOLDH»2/2. 

1680=  H I NEU=C2« ( < VRHO- VNEH 1 /VSVNEH1 4*2/2 . 

1690=  P(I)  = (PUU(Z.*VOLD-CAHMA4VDOT)-2.»GANMA*(HINEW-H10LD)4 

1700=  12.*(VNEW4HPNEH-VOLD»HPOLOI*2.4CAMMA»0(I)*VDOT1/ 

1710=  2(2.4VNEW4GAHHA4VD0T) 

1720=  GO  TO  40 

1730=  38  CONTINUE 

1740=C COMPUTE  PRESSURE  USING  VARIABLE  GAMMA  

1750s  S*2./(VC(I)-1.) 

1760=  Pls(P(IH(S4VOLD-VDOT)-2.4Q(I)»VDOTI/fS»VNEW4VDOT) 

1770=  EI=EIA(PliDtI)»RHOI 

1730=  VG ( I ) = VGAMMA (E I « D ( I ) » RHO) 

1790=  Ss2./(VG(I)*1.1 

1800=  P(I)=(P(l)i(S»VOLO-VDOT)-2.*8<Il4VDOT)/(S*VNEW4VDOT) 

1810=  40  CONTINUE 

1820=C LOCATE  NEW  SHOCK  POSITION --- 

1830*  NP5=N145 

1840=  NH5*N1*5 

1850*  DO  52  I=NM5iNP5 

1860*  IF(OMAI.CT.QU))CO  TO  52 

1870*  QMAX*Q(I) 

1880*  N1=I 

1890*  52  CONTINUE 

1900=C ADVANCE  TIME 

1910*  TIME*TIME4DT 

1920 *C HTDROHATIC  TINE  STEP  CONTROL 

1930*  IF (NSAVE-Nt .EO.01RETURN 

1940*  1F(IC0UHT.LE.51TFRACS.5»T"RAC 

1950*  IF  UCOUKT .GT. 101 TFRAC*TFRAC+. I 

1960*  ICOUNT*0 

1971*  RETURN 

1981*  END 


l 


1990* 

,4 

\ 

2000* 

SUBROUTINE  SQUEEZE (CA *L. CAMMA.N1 »DI » RHO.P.D. X ,Q(U»CZfDT»CK) 

h 

2010* 

DIMENSION  PU0l)tDli0l). 1(101)  >0(101). U(  101)  »CK(7) 

i , 

2020* 

CA*4 

>i 

2030sC- 

READ  PROBLEM  GEOMETRY ,CA«HA»DENS1TY  tSOUND  SPEED. 

i 

2040*0 

RADIUS  OF  METEL.  tND  HASS  OF  HIGH  EXPLOSIVE  

; i 

2050* 

READ*. L.CAHMA.RHO.C. RADIUS. PRESS 

2060s 

HRITEtiL.CAMMA.RHQ.C. RADIUS. PRESS 

j 1 

2070=C- 

COMPUTE  INITIAL  CONDITIONS  FOR  PROBLEM  

2080* 

DX*RADIUS/I00. 

2090* 

C2=C**2 

i i 

2100* 

PRINT*."DXs".OX."  AND  C2*".C2 

2110* 

PNORM=1.014E6 

2120*C- 

ASSIGN  SCALED  VALUES  TO  PROBLEM  ARRAYS  

i 

2130* 

DO  10  1*1.101 

2140* 

U(I)*Q(I>*0 

' i 

2150* 

P(II=PNORH 

i1 

2160* 

Dtll'RHO 

1 

2170* 

IF (I .EQ. 1) X (I) =0 

■ If 

2180* 

IF(I.CT.miI)*I(H)*DI 

|j 

2190* 

10  CONTINUE 

- ! 1 

ZZ00=C INITIALIZE  LOCAL  PARAMETERS 

2210=  HNaNl*100 

2220=  OPOS*X(Nl)  i 

2230s  TlHIT=CtCLEsKFIFT^=§  f 

2240s  P(100)*PRESS 

2250-  RETURN  . ; ! 

2260*C  j j 

2270s  ENTRY  CRUSH  | 

2280=C INPUT  A CARO  HERE  FOR  CONST  APPLIED  PRESSURE ; ' 

2290SC CARD  SHOULD  READ  (PUOdhPRESS)  I f 

2300s  IF (P ( 100) .LT.PHORH) P 1 100} =PNORH  ! 

2310*0 RESET  SEMETRIC  U.  C.  AND  ADVANCE  TIME i 

2320s  XU)=U(l)=0.  -j 

2330*C TIME  SHOCK  PROGRESS  FOR  50  ITERATIONS , 

2340*  KFIFTf*KFIFTT+l  - I 

2350s  TINIT*T1NIT+DT  \ \ 

2360*0 COMPUTE  OVERPRESSURE  AND  OVERDENSITY  1 

2370=C PLACE  COMPUTED  VALUE ^ IN  CHECK  ICK)  ARRAY j j 

2380*  IP=N1  i 

2390*  CK(IJ*CTCLE*CTCI.E*1.  j 


2400*  CKI2)*T1HE*TIME*DT 

2410*  CKI3)*DT 

2420*  CK(4)*PflS»X(IP) 

2430*  PHAX‘DMAMMAI‘0. 

2440*  NH5*Nl-5 

2450*  NP5*Nl+5 

2460*  DO  20  l*NN5iNP5 

2470*  IF (VMAI .LT .U ( It  I VMAI*U( II 

2480*  lFIP(I) .GT.PMAl)PNAI*P(II 

2490*  IF(D(I).GT.DHAI)DI1AI*D(Ii 

2500*  20  CONTINUE 
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2510*  CK (5) =OP=PHAI“PNORM 

2528*  CK {&) =OD=DMftX-RHO 

2538 =C COMPUTE  NEU  MATERIAL  VELOCITY  

2548*  CK(7)*VMAI 

2558*  IF (KFIFTY.LT.SflICO  TO  48 

2568*  OPOS*PO$ 

2578*  TIHlT=KFIFTT=0 

2580*C WRITE  CHECK  VALUES  EVER!  58  ITERATIONS 

2598*  URITE(8t2000l(CK(l)iI*lt>) 

2688*  48  COIITIHUE 

2618*  110  FORMAT!/. 101. "CORE  IS  SOHO  WITH  ".C18.2."  CHS  HE."./) 

2628*  2800  FORMAT ( IX ."CYCLE*" . F5.0 T!SEC)*''.IPE9.2."  DT(SEC)*" 

2630=  1.E7.0."  POS(CMI*"»E9.2i"  0P!D/CH**2)*".E9.2»"  0D(G/CM»*3)*« 

2640*  2.E9.2."  MVEL(CM/SEC)=".E9.2) 

2650=  RETURN 

2660*  END 

2670= 

2660*  SUBROUTINE  BLAST (CA.L.GAMMA .Ml .M.DIiRHO.P.D.X.Q.U.DT.CK) 

2690*  DIMENSION  PI101)  .01181)  .1(101)  .0(161)  iU(101)  .CK(7) 

2780=  SCALE ( INQRM i X D ) * INORM*  T KT** ( 1 . /3 . ) * ID 

2710=  CA=4 

2720*C READ  YIELD.SCflLIHC  PARAMETERS. AND  PROBLEM  GEOMETRY 

2730=  READMKT.SP.SD.ST 

2740*  L*3 

2750=  WRITE* .YKT.SP.SD.STjL 

Z760=C COMPUTE  INITIAL  CONDITIONS  USING  SCALING  LANS ' 

2770*  PNORM=1.014E6*SP 

2780*  CK(2)*TIME=SCALE(3.68E-5tST) 

2790*  PBLAST*1.5E10»SP 

2800*  DI=SCALE(50.iSDI 

2810*  RH0=1.293E*3*$D 

2820*  Nl=9 

2830*  0P0S=I(N1) 

2840*  PRINT* ."PNORH*" . PNORM i " PBLfiST="»PBLAST 

2850*  PRINT*."  TIME*".TIME."  DI*".DI."  RHO*"»RHO."  Nl«".Nl 

2860*C ASSIGN  SCALED  VALUES  TO  PROBLEM  ARRAYS  

2870*  DO  18  1*1.111 

2880*  QU)=U(I)=0 

2890*  D(1)=RH0 

2900*  IFtl.EQ. 1)1111=0 

2910*  IF(LGT.l)im*I(I-l)*OI 
2920*  IF(I.LE.Nl)P(I)=PBLAST 

2930*  1F(I.CT.HI)PU)*PN0RH 

2940*  10  CONTINUE 

2950*C- INITIALIZE  LOCAL  PARAMETERS  

2960*  N*Nl+5 

2970*  T1H1T*CTCLE*KFIFT!*0 

2980*  RETURN 

.2990* 
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3BM*  ENTRY  EXPAND 

3010=0 RESET  SEMETRIC  B.  C.  AND  ADVANCE  TIKE 

3020*  Xll)=U(l>=0. 

3030*C LOCATE  NEH  SHOCK  POSITION 

3040*  IP=N1 

3050*C TIKE  SHOCK  PROGRESS  FOR  50  ITERATIONS 

3040*  KFIFTT*KFIFn+l 

3070*  TIH!T*TINITH)T 

3080*C COMPUTE  OVERPRESSURE  AND  OVERDENSITT 

3090=  PHAI=DMAX=VMAX=6. 

3100*  NM5=Nl-5 

3110*  NP5=Nl+5 

3120*  DO  20  I=NK5iNP5 

3130*  IF(UII) .GT, VMAX)VMAX*U( I ) 

3140*  IFIP(I)  .CT.PKA1)PMAI(*P(I) 

3150*  IFID(I) .GT.DHAX)DNAI*D(1J 

3160*  20  CONTINUE 
3170*  CK(5)=0P*PMA!i-FN0RM 

3180*  CK(6)*0D*DNA)t-RH0 

3190=C PLACE  COMPUTED  VALUES  IN  CHECK  iCK)  ARRAY 

3200*  CK(1)*CYCLE*CTCLE+1. 

3210*  CK(2)vTIHE*TIME4DT 

3220*  CK(3I=DT 

3230*  CK(4)*P0S*X(IP) 

3240*C COMPUTE  NEW  SHOCK  MATERIAL  VELOCITY 

3250*  CK(7)*VMA2 

3260*  IFIKFIFTY.LT. 50ICO  TO  30 

3270*  OPOS*POS 

3280*  TIMIT*KFIFTY*0 

3290*C WRITE  CHECK  VALUES  EVERY  50  ITERATIONS 

3300*  NRITE(8i2000MCK(I)tI*li7) 

3311*  URITE(2>100HCK(K)'K*1'7) 

3321s  3f  K-NU5 

3330*  100  FORMAT (21i"CYCLE*"iF6.0 i"  TIME*"tlPE9.3."  DT*"iE9.3. 

3340*  l"  S POS  *"iE9,3i"  0 PRES«”iE9.3. " 0 DEN*"E9,3»*'  S VEL*" 

3350*  2iE9.3) 

3360*  2000  FORMATUXi"CTCLE*"tF5.0i"  T1SEC)*MPE9.2»H  DTISEC)*" 

3370*  1iE7.0i"  P0S(CH)*"iE9.Zi"  OPIO/CMm2»*".E9.Z«"  0DIG/CH«i3)*" 

3380*  2.E9.2."  HVEL(CM/SEC)*M9.2> 

3390*  RETURN 

3400*  ENO 

341f*C 

3420*  SUBROUTINE  REZONEIXiUiPi D>Q>DX iCAHMAi ICYCLE»TIHE*N1 >RHO) 

3430*  DIMENSION  X (If l)  fU(191) »P(101) tDUf  11  >PUB<  101  *3) »OU01 ) 

3441*  Nl*Nl/2 

3451*  WRITE(2»100)TlHEt ICICLE 

3460*  DO  10  J*3»97i2 

ii70*  I*W*l)/2 
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3480*C COMPUTE  GAMMAS 

3490=  EIJl=EIA(P(JltO(J).RHOI 

3500=  EIJ2=EIA(P(J+1) iDIJ+t) iRHQ) 

3510=  AGAMA=VCAMMA(Ei JliD(J) tRHO) 

3520*  BGAMA'VCAMHA (El J2 1 0 ( J+ 1) i RHOI 

3530=C COMPUTE  HASS  AMD  NEW  DENSITT  

3540=  XMJ1=D  ( JI*<X  ( J+H  **3-K  < Jl  **3) 

3550=  XH.I2=DlJ+l)*(X(J+2)*«3-X<J+lMt*31 

3560«  PUD(Ii3)=(XHJl+XHJ2)/(X(J+2U*3-XtJ)*«3> 

3570o  XMI1=PUD  (I»3)*  (X  (J+Z)*«3-X<J)-»*3) 

3580=C COMPUTE  OLD  CELL  CENTERED  VELOCITY 

3590*  UJl=(U(J)+U{J+l))/2. 

34000  UJ2= IU ( J+i) +U(J+2) ) /2. 

3616*C COMPUTE  OLD  ENERGIES  (KENETICt  INTERNAL)  TOTAL)  — 

3620*  EKJl=UJU*2/2. 

3630*  EKJ2oUJ2H2/2. 

3640*  EIJi*P(J)/({ACAHA-l.MD(JI) 

3650*  E1JZ=P(J+1)/((EGAHA-1.MD(J+1)) 

3660=  ETJloEKJHElJl 

3670o  ETJ2=EKJ2+EIJ2 

’ 3680*C COMPUTE  NEW  TOTAL  ENERGT 

3690*  ETI1=(XMJUETJ14XMJ2*ETJ2)/XMI  1 

3700*C COMPUTE  MOMENTUM 

37t0o  XHONJt=XMJHUJl 

3720*  XH0MJ2*  XM J2*U J2 

3730*  PUD(  1 .2)  = (XM0NJ1+XN0MJ2)  /INI  l 

3740*  XHOHI 1=XKI!*PUD ( I »21 

3750=C COMPUTE  NEW  KENET1C  AND  INTERNAL  ENERGT 

3760o  EKI|oPUD(IiZM*Z/2. 

3770o  EIll'ETIl-EKIl 

3780«C COMPUTE  NEW  PRESSURE 

3790o  GAIU1A'VGAMMAlElIliPUB(l»3)  iRHO) 

3800*  PUDU.l)*Eni»(CAMMA-U*PUDIIi3l 

3810*C BTPASS  LENGTH T REZONE  PRINTOUT 

3820*  IF (Ml .GY . 10) GO  TO  10 

3830*  tlR!TE(2i200)X(J)  tUJliPlJ)  iD(J)«ftlJliXMOHJlfETJl*EKJi>ElJl 

3840o  WRnE(2i200)X(J+l)>U<J2iPU+l)>D(4dl»XHJ2iXMOHJZiETJZ>EKJ2iElJ2 

3850o  URITE(2i250)X(JhUJIiPUCUil)iPUD(I)3IiIHIhXMOHIiiETUiEKIliEIIl 

3860*  10  CONTINUE 

3870*C ASSIGN  VELOCITT  FOR  CENTER  AND  RESET  DE 

3880*  PUD(I)2)*0 

3890*  DX=2.*DX 

— STORE  NEW  VALUES  IN  PROGRAM  ARRATS  

DO  20  I*li49 
J*2»I-1 
X(I)oXIJ) 

PU)*PUO(l»l) 

Dm*PUDII.3» 

IFd.EQ.DCOTO  20 


O 


v 

3910* 

3920* 

3930* 

3940« 

3950* 

3960* 
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3970*C COMPUTE  NEW  BOUNDARY  VELOCITIES 

3988s  U(I)*(PUD(I-ltCI+PUD(Ii2) 1/2. 

3990s  20  CONTINUE 

4000SC DEFINE  VALUES  FOR  UPPER  HALF  OF  ARRAYS 

4010s  DO  25  I*50i 101 

4020s  MbUMHDI 

4030s  U(I1*Q(I)=0 

4040s  P(1)*P(1G1) 

4050s  D(I)=D(101) 

4060s  25  CONTINUE 

4070*0 COMPUTE  NEU  VISCOSITY  - 

4080*  DO  30  1=1 »50 

4090*  Q<I)=8.*D(I)*(UIH1)-UII))m2 

4100*  IF(UU+1)-U(I).CE.0IQ(1)=0 

4110*  30  CONTINUE 

4120*C SET  BOUNDARY  CONDITIONS 

4130*  U<1)=X(1>=0 

4140=  PU1=P<2) 

4150*  Q(1)*Q<2) 

4160*  D(1)=D(2) 

4170*  100  FORMAT  I 1H1iT20i"KE2ONE  TIME  * "iC9.3."  ITERATION  NUMBER  M7i 
4180=  l////»T4i"P0SITI0N"iT18i"VEL0CITT,'iT32l"PRESSURE"iT46i "DENSITY" 

4190=  ZT&2* "MASS" t T74 » "MOMENTUM" » T80 1 -'T  EMERGY"iT102t,,K  ENERGY"! 

4200*  3Tll6i"I  ENERGY"  iM7i"CM"iT19i"CM/S£C"iT33i"D/CNh3"» 

4210*  4T46."G/CHH3",T63i"GMl,iT76i"C-SEC"»T91."ERC/G"iT105»"ERG/GM. 

4220=  5TI19i"ERC/C"i///) 

4230s  200  FORMAT (3Xi9(G9. 3.5D) 

4240*  250  F0RHATI51t.9(G9.3.5*)iYI 
4250=  RETURN 

4260*  END 

4270* 

4280=  FUNCTION  EIAIPiIi.RHO) 

4290*CC GUESS  GAMMA  MINUS  ONE 

4300=  GM1*.2B 

4310=  DO  30  I=lil0 

4320=0 COMPUTE  AN  INTERNAL  ENERGY  

4330*  tI=P/ICNHDl 

4340*C CONFUTE  A NEH  GAMMA  MINUS  ONE  

4350*  CHI *VGANHA IE I »D iRHOl - 1 . 

4360*C COMPUTE  A PRESSURE 

4370*  PZ*GH1*D*E1 

4380*C COMPARE  PRESSURES - 

4390*  IF(ABS(P-P2).LE,.001iPIGO  TO  40 

4400*C GAMMA  IS  CORRECT  WHEN  THE  PRESSURES  AGREE 

4410*  30  CONTINUE 

4420*  PRINT*»"l*M 

4430*  40  EIA*EI 

4440*  RETURN 

4450s  END 

1460* 


O C. 


{ 


4478=  FUNCTION  VCAHMA(EI,D,RHO) 

4480-C HUFF  ROUGH  VARIABLE  CAMI1A  COMPUTATION 

4490=  EIE=EI/1.E10 

4500=  IF (EIE.CT. . L315IG0  TO  10 

4510=  CN0NE=.398l 

4520=  CO  TO  3fl 

4530=  10  CONTINUE 

4540=  IFlEIE.CT.l.lGO  TO  20 

4550=  CNONE=-.0399MEIE-.1315H»2+.3981 

4560=  CO  TO  30 

4570=  20  CONTINUE 

4580=  GH3NE=. It*. 624/(2. +EIE) 

. 4590=C ADJUSTMENT  FOR  DENSITY  

4600=  RHQZ=D/RHO 

4610=  EL=ALOC10(EIE) 

4620=  ALPIIA=.0577»EL-.C218iELH2+.0035»EL4f3+.0002*Em4 

4630=  CHQNE=GHONE*RHOZ**ALPHA 

4640=  30  VGAMNA=CHONE+l« 

4650=  RETURN 

4660=  END 

4670= 

4680=  SUBROUTINE  COPY  IX *U»P.DiQ»FLAG»TIME*TAELE»CKl 

4690=  DIMENSION  I ( 1 01 ) >U( 101 ) iP( 101 ) >D ( 101 ) iQ ( 101 ) : IABLE ( 10) 
4700=  liCK(7) 

4710=  INTEGER  FLAG 

4720=C RESET  CALLING  PARAMETERS 

4730=  FLAG=FLAG-1 

4740=  DO  10  1=1.9 

4750=  10  T ABLE ( I ) =T ABLE ( I + 1 ) 

4760=  1F(FLAG.LT.0ITABLE(1 1=1000. 

4770=C WRITE  INFORMATION  ON  LOCAL  FILE  (FILES)  

4780=  ENTRT  ICOPT 

4790=  NR!TE(7.10{i0)TIHE.  (XU),  U(I).P(I).D(I).QU).  1=1.101) 

4H0=  NR1TE(7.2000I(CK(K).K=I.7) 

4810=  WRITE(2i1010)T1HE 

4820=  CALL  PLOTSHO(.P.100I 

4830=  NRITE12.2000) (CK1K).K=1.7) 

4840=  URITE  <£• 1020) TIME 

4850=  CALL  PLOTSMil.D. 100) 

4860=  WRITE (2,2000) (CK(K) iK=1.7) 

4870=  WR1TE(2,1001)TIME 

4880*  CALL  PLGTSMU.U,  100) 

4890=  NR1TEI2, 20001 (CK(K), K=l, 7) 

4900=  WRITE(2,1030ITIME 

4910=  CALL  PLGTSHd, Q.100! 

4920=  URITE(2,2000) (CK(K) , K* 1 ,7) 

O 
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1030  FORMAT (IX t"TlKE  IS;,.F9.5."SEC".//.T7,"P0S!T10N".T22 
4940s  li"VEL0Cm"»T37.' "PRESSURE".  T52i  "DENSITT"»Tt7»,,VISCQSITT" 

4950=  2»/»T9."!C[II"»T22."(CH/SEC)".T37."!D/CHh2)".T51."(G/CMh3)" 

4960=  3»T67t"(D/CMH2)".//.10l  (5(5XilPE10.3l  ./)) 

4970=  1001  FORMAT ( 1H1  > T40 1 "VELOCITY !CM/S£C)  VS  POSITION(CM)  (IT  TIME  = " 
4980=  1 i 1PE10 .3 i " SEC") 

4990=  1010  FORMAT!  Ml  it40."PRESSURE!D/CHh2)  VS  POSITIOM(CM)  AT  TIME  = " 
5000=  1i1PE10.3i"  SEC") 

5010=  1020  FQRMATllHliT40."DENSlTT !G/CH»*3)  VS  POSITIOH(CM)  AT  TIME  = "» 
5020=  11PE10.3."  SEC") 

5030=  1030  FORMAT! llll tT40t"VISCOSlTl  (D/CM*«2>  VS  POSITIGN(CM)  AT  TIKE  = " 
5043=  1.1PE10.3."  SEC") 

5050=  2000  FORMAT !l)(i"CTCLE="iF5,0i"  T(SEC)=".IPE9.2."  BT(SEC)=" 

5060=  liE7.0»"  P0S(CMI=">E9.2i"  0P(D/CHh2)=".E 


5060=  1.E7.0."  P0S(CMI="iE9.2»"  0PiD/CH«2)=".E9.2."  OD(C/CMm3)=" 

5070=  2.E9.2."  MVEL(CH/SEC)="iE9.2) 

5080=  RETURN 

5090=  END 

5100= 

5110=  SUBROUTINE  PLOTS«(X»Y ,H) 

5120=  DIMENSION  1(100)  .Y(100) , III  100.40)  .HLI40) 

5130=  NP1=N+1 

5140=  DATA  SLASH . DOT . ST AR i OH . BLANK/ 1 H/.1H . » 1 H* . 1HQ • " "/ 

5150=  XMAX=X (N+l I 

5160=  TMAI=TKINr0 

5170=  DO  10  1=1 iN 

5180=  IF(T(I) .CT.TMAX) YMAX=T ( I) 

5190=  IF(T(H.LT.TMIN)THIN=TII) 

5200=  10  CONTINUE 

5210=  TOVER=THAI-T!100I 

5220=  DO  20  1=1.100 

5230=  DO  20  J=1 »40 

5240=  U(I.J)=BLA)IK 

5250=  IF(I.E0.1)W(I.J)=SLA$H 

5260=  IF(J.EO.20)H(I,J)=DOT 

5270=  H(10.20)=N(20t20)=H(30i20)=H(40i20!=H(50)20)=H(99.20)="0" 

5280=  H(60i20)=H(70.20)=H(80i20)=H(90i201=NI100»20)="0" 

5290=  H(V.20)=U(98iZ0)="l" 

5300=  H(19i20)*"2" 

5310=  H(29.20)="3" 

5320=  11(39. 20)="4" 

5330=  11(49. 20)="5" 

5340=  11(59. 20)="6" 

5350=  H<69.20)="7" 

5360=  «<79,20)="8" 

5370*  H<89.20)="9" 

5380=  20  CONTINUE 


5390=  SCALE-TNAX/20. 

5400=  STORE=-tl1IM/20. 

5410=  IF (SCALE. LT. STORE) SCALE=STORE 

5420=  DX-XKAX/5. 

5430=  SCALE=$CALE+.1»SCALE 

5440=  G«AX=20.«SCALE 

5450=  DO  30  1=1 .21 

5460=  X1K1= ( I- 1 ) 

5470=  ST0RE=X1H1«SCALE 

5480=  K=21-I 

5490=  L=19+I 

5500=  HL(L)=-STORE 

5510=  WL (K) =STORE 

5520=  30  CONTINUE 

5530=  DO  40  1=1.26.5 

5540=  J=4»  (1-1) 

5550=  n=(I-l)/5 

5560=  UL(1)=X1»DX 

5570=  40  CONTINUE 

5580=  DXX=XMAX/100. 

5590=  DO  50  1=1. N 

5603=  DO  44  K=1.N 

5610=  TK=K 

5620=  IF(X(I).GT.XK*DKX-DIX/2..AHD.X(I).LE.XK*DXX+DXX/2.)G0  TO  45 

5630=  44  CONTINUE 

5640=  45  CONTINUE 

5650=  DO  50  J=1.40 

5660=  XJH1=J-1 

5670=  XJ=J 

5680=  IF  ( (GKAX-XJMHSCALE.CT.  T < I ) ) .AMD.  (CHAX-X  J«SCALE 

5690=  i.lE,T(I)))U(K.JI=OH 

5700=  50  CONTINUE 

5710=  UR1TE(2» 1001XHAX.XMAX 

5720=  K=5 

5730=  WRITE (2.200) GhAI 

5740=  DO  60  J= 1.40 

5750=  IF < J.EQ.K) GO  TO  55 

5760=  VRITE(2. 400101(1. J).I=1»NI 

5770=  CO  TO  60 

5780=  55  K=K+5 

5790=  WRITE (2i200)UL< J) . (H(J. J1 .I=i.N) 

5800=  60  CONTINUE 

5810=  WRITE (2.360) (WL(I). 1=1.26. 5) 

5820=  100  FORMAT (y // . T30. "THIS  PLOT  HAS  I«A1=  MPE10.3.".THAI(=  " 

5830=  l.E10.3.MHIN=  ".E10.3."  AND  t OVER*  ".E10,3.//y) 

5840=  200  FORMAT  I5X.G10.3.2M06A1) 

5850=  300  FORMAT (/// .T14.6I 1PE10. 3. 10X)) 

5860=  400  FORMAT (17X. 100AL) 

5870*  RETURN 

5880=  END 
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Appendix  F 

User's  Guide  for  HUFF 

This  appendix  describes  the  data  card  input,  the  scope  control 
cards,  and  selecting  the  output  when  operating  program  HUFF.  The  data 
card  input  will  be  presented  first  followed  by  a discussion  of  the  scope 
control  cards.  The  last  section  will  discuss  the  types  of  output  avail- 
able to  the  user. 

Data  Cards  for  HUFF 

HUFF  is  programmed  to  solve  two  types  of  problems.  The  first  type 
of  problem  is  the  nuclear  blast  in  air  and  the  second  type  involves  the 
compression  of  solid  materials  with  a constant  applied  pressure.  Both 
problems  require  on  input  of  three  data  cards.  The  first  two  data  cards 
establish  the  parameters  used  in  operating  the  programs  and  will  be  dis- 
cussed first.  The  last  card  is  the  problem  defining  card  and  has  a for- 
mat based  on  the  type  of  problem.  All  the  data  is  read  by  HUFF  with  a 
list  directed  format  allowing  the  user  maximum  flexibility  in  defining 
the  input  numbers . 

Parameters . The  two  input  cards  defining  the  program  parameters 
are  as  follows! 

1.  Card  1 - FLAG  (2),  FLAG  (3),  FLAG  (4),  TSTOP . This  card  con- 
tains the  input  flags  defined  as  follows; 

, FLAG  (2)  is  an  integer  representing  -the  number  of  result 
timeB  listed  on  card  2. 

FLAG  (3)  is  an  integer  representing  the  maximum  number  of 
time  step  iterations . The  program  will  Btop  if  it  reaches 


90 


this  many  iterations . 


1 


I 


I 


i 

£ 
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FLAG  (4)  is  an  integer  representing  the  type  of  problem  to  be 
solved 

1 for  a problem  in  air 

2 for  a problem  in  a solid 

TSTOP  is  the  maximum  problem  time  allowed.  The  program  will 
stop  if  it  reaches  this  time.  Results  for  this  time  will  be 
plotted  in  addition  to  those  requested  on  Card  2. 

2.  Card  2 - TABLE  (1),  TABLE  (2),  . ..,  TABLE  (FLAG (2)).  This  card 
lists  the  requested  result  times.  A maximum  of  10  times  can  be 
requested  due  to  the  dimension  of  the  array  TABLE. 

Problem  Defining  Card.  This  is  the  third  and  last  input  card  needed 
to  define  the  problem.  The  blast  problem  will  be  discussed  first  followed 
by  the  compression  problem. 

1.  Card  3 - YKT,  SP,  SD,  ST  is  the  problem  defining  card  for  the 
blast  problem.  The  variables  are  defined  as  follows; 

YKT  - Yield  of  the  blast  in  kilotons . 

SP,  SD,  ST  - Altitude  scaling  factors  from  Table  I on  page  15 
of  the  text. 

2.  Card  3 - L,  GAMMA,  RHO,  C,  RADIUS,  PRESS  are  the  problem  de- 
fining variables  for  the  compression  problem.  These  variables 
are  defined  as  follows : 

L - Problem  geometry  factor 

1 - Rectangular  geometry 

2 - Cylindrical  geometry 

3 - Spherical  geometry 
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GAMMA  - Gruneisen  ratio  for  the  solid 


RHO  - Normal  density  of  a solid 

C - Normal  density  sound  speed  in  the  solid 

RADIUS  - Uncompressed  radius  of  the  solid 

PRESS  - Pressure  applied  in  the  problem  by  the 
outer  cell 

The  pressure  in  the  compression  problem  can  be  applied  throughout 
the  problem  solution  by  including  the  card  (P(100)=PRESS)  after  ENTRY 
CRUSH  in  subroutine  SQUEEZE.  If  this  fortran  card  is  not  included, 
the  pressure  will  be  applied  only  in  the  initial  conditions  corresponding 
to  an  impulse  pressure . A comment  card  is  included  at  the  proper  loca- 
tion in  subroutine  SQUEEZE  to  identify  this  option. 


Scope  Control  Cards  for  HUFF 

This  section  presents  the  scope  control  cards  necessary  to  execute 
the  HUFF  program  at  AFIT.  The  sequence  of  control  cards  used  will  depend 
on  what  the  user  plans  to  do  with  the  program  and  its  results  , A number 
of  options  will  be  discussed  and  explained  to  provide  flexibility  for 
the  user. 

Permanent  File . To  place  HUFF  on  permanent  file  and  execute  a pro- 
gram simultaneously,  the  following  control  cards  are  necessary: 

job,T200 .problem 
REQUEST , HUFF , *PF . 

COPYBR , INPUT , HUFF . 

REWIND ,HUFF . 

CATALOG ,HUFF , RP-999 , 

REWIND, HUFF. 

FTN , I-HUFF , 

LGO. 

7/8/9 
HUFF  DECK 
7/8/9 

DATA  Cards  for  HUFF 
7/8/9 

6/7 /8/9  (Orange  Card) 
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After  HUFF  1b  on  permanent  file,  the  following  control  cardB  can 


be  ueed; 


( 


job , T200 .problem 
ATTACH ,HUFF . 

FTN.I-HUFF. 

LGO. 

7/8/9 

DATA  Cards  for  HUFF 
7/8/9 

6/V/8/9  (Orange  Card) 

Single  Operation . The  following  cards  use  HUFF  to  3olve  a single 

problem  without  using  a permanent  file : 

job ,T200, problem 
FTN. 

LGO. 

7/8/9 

HUFF  PROGRAM 
7/8/9 

DATA  for  HUFF 
7/8/9 

6/7/B/9  (Orange  Card) 

Listing  Local  Files . HUFF  writes  all  the  plotted  information  and 
summary  information  on  local  files  called  FILES  and  REVIEW  respectively. 
The  following  cards  will  route  a listing  of  the  local  files  to  the 
printer : 


job,T200 .problem 
FTN. 

LGO. 

REWIND .FILES. 

ROUTE  .FILES  ,TID*BB  ,DC«PR,FID»  job  . 

REWIND, REVIEW. 

ROUTE, REVIEW, TID*BB,DC=PR, FID- job. 

7/8/9 

HUFF  PROGRAM 
7/8/9 
DATA 
7/8/9 

6/7/8/9  (Orange  Card) 

Save  Local  Files.  To  save  local  files  on  a permanent  file,  the 
following  cards  are  necessary « 
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job ,T200 .problem 

Request , If n , *PF . 

FTN. 

LGO. 

REWIND  ,lfn . 

CATALOG, If n,RP-9 99. 

7/8/9 

huff  program 

7/8/9 

DATA  for  HUFF 

7/8/9 

6/7/8/9  (Orange  Card) 

There  are  many  more  combinations  available  in  using  scope  cards  to 
run  a program.  The  above  listings  should  provide  sufficient  examples 
to  operate  the  program. 

HUFF  Output . A large  amount  of  useful  information  1b  generated 
during  a problem  solution.  This  information  is  easily  analyzed  if  it 
is  presented  in  graphic  form.  However,  if  a more  detailed  analysis  of 
the  results  is  required,  a tabular  output  may  be  preferred.  For  this 
reason,  a variety  of  outputs  are  generated  by  HUFF.  A tabular  form  of 
all  the  results  is  placed  on  a local  file  called  FILES , a summary  of 
the  parameters  of  the  problem  is  placed  on  a local  file  called  REVIEW 
and  all  the  information  is  plotted  on  the  output  file  by  subroutine 
PLOTSM. 

Plot  routines  using  DISSPLA  or  CALCOMP  can  also  be  used  with  HUFF 
by  replacing  subroutine  PLOTSM  with  the  desired  plot  routine.  If  FILES 
and  REVIEW  are  placed  on  a permanent  file,  selected  data  can  be  read 
from  these  files  and  plotted,  after  examing  the  results  displayed  on 
the  output  file,  without  reexecuting  the  entire  program. 
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HUFF  is  a one-dimensional  Lagrangian  hydrodynamics  computer  code  developed] 
from  the  basic  principles  of  mass,  momentum,  and  Energy  conservation  for 
strong  shock  propagation  in  a solid  or  gas.  Two  equations  of  state  are  used 
- the  adiabatic  ideal  gas  law  with  a variable  gamma  and  the  Gruneieen  solid 
equation  of  state  with  a cone tent  Gruneieen  ratio.  The  Richtmyer  and  Morton 
difference  equations  for  strong  shocks  are  used  on  a spatial  mesh  composed 
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a usefulness  and  limitations  of  the  code  and  sIbo  serve  as  sample  problems. 
The  results  of  a one  kiloton  nuclear  explosion  are  compared  to  the  Nuclear 
Blast  standard  1KT . The  results  wore  within  13  percent  for  shock  over- 
pressure and  overdensity,  5 percent  for  shock  material  velocity,  and  2 
4 ircent  for  shock  position  over  a range  of  20  meters  to  2 kilometers  front 
the  burst  point.  The  larger  deviations  occurred  at  early  times  being 
attributed  to  an  absence  of  radiation  transport  calculations  in  the  code . 

The  second  problem,  a megabar  compression  of  uranium,  shows  agreement  within 
two  percent  for  all  parameters  (peak  shock  pressure,  density,  material  velo- 
city and  shock  velocity)  when  compared  with  the  Rankine-Hugoniot  compression 
curves.  The  equation  of  state  for  a solid  was  limited  to  calculations  below 
100  megabars  due  to  its  simplicity  and  constant  value  for  the  Gruneieen 
ratio.  A complete  users  guide  and  program  listing  are  also  provided. 
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